diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index c241153e76..fcabd989e0 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -162,6 +162,7 @@ global.o: cmfd_header.o global.o: constants.o global.o: dict_header.o global.o: geometry_header.o +global.o: hdf5_interface.o global.o: material_header.o global.o: mesh_header.o global.o: particle_header.o @@ -171,17 +172,16 @@ global.o: source_header.o global.o: tally_header.o global.o: timer_header.o -hdf5_interface.o: ace_header.o -hdf5_interface.o: bank_header.o -hdf5_interface.o: constants.o -hdf5_interface.o: endf.o -hdf5_interface.o: error.o -hdf5_interface.o: geometry_header.o -hdf5_interface.o: global.o -hdf5_interface.o: material_header.o -hdf5_interface.o: output.o -hdf5_interface.o: string.o -hdf5_interface.o: tally_header.o +hdf5_summary.o: ace_header.o +hdf5_summary.o: constants.o +hdf5_summary.o: endf.o +hdf5_summary.o: geometry_header.o +hdf5_summary.o: global.o +hdf5_summary.o: hdf5_interface.o +hdf5_summary.o: material_header.o +hdf5_summary.o: output.o +hdf5_summary.o: string.o +hdf5_summary.o: tally_header.o initialize.o: ace.o initialize.o: bank_header.o @@ -262,6 +262,9 @@ output.o: plot_header.o output.o: string.o output.o: tally_header.o +output_interface.o: hdf5_interface.o +output_interface.o: mpiio_interface.o + particle_header.o: constants.o particle_restart.o: bank_header.o @@ -336,8 +339,8 @@ source.o: string.o state_point.o: error.o state_point.o: global.o -state_point.o: math.o state_point.o: output.o +state_point.o: output_interface.o state_point.o: string.o state_point.o: tally_header.o diff --git a/src/OBJECTS b/src/OBJECTS index c9612e6af0..b079906751 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -29,6 +29,7 @@ geometry.o \ geometry_header.o \ global.o \ hdf5_interface.o \ +hdf5_summary.o \ initialize.o \ interpolation.o \ input_xml.o \ @@ -38,7 +39,9 @@ material_header.o \ math.o \ mesh_header.o \ mesh.o \ +mpiio_interface.o \ output.o \ +output_interface.o \ particle_header.o \ particle_restart.o \ particle_restart_write.o \ diff --git a/src/constants.F90 b/src/constants.F90 index e9688e761d..4442d05ca1 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -11,8 +11,13 @@ module constants integer, parameter :: VERSION_RELEASE = 1 ! Revision numbers for binary files - integer, parameter :: REVISION_SOURCE = 1 - integer, parameter :: REVISION_STATEPOINT = 8 + integer, parameter :: REVISION_STATEPOINT = 9 + integer, parameter :: REVISION_PARTICLE_RESTART = 1 + + ! Binary file types + integer, parameter :: & + FILETYPE_STATEPOINT = -1, & + FILETYPE_PARTICLE_RESTART = -2 ! ============================================================================ ! ADJUSTABLE PARAMETERS @@ -359,10 +364,9 @@ module constants integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file integer, parameter :: UNIT_XS = 14 ! unit # for writing xs summary file - integer, parameter :: UNIT_SOURCE = 15 ! unit # for writing source file - integer, parameter :: UNIT_STATE = 16 ! unit # for writing state point - integer, parameter :: CMFD_BALANCE = 17 ! unit # for writing cmfd balance file - integer, parameter :: UNIT_PARTICLE = 18 ! unit # for writing particle restart + integer, parameter :: CMFD_BALANCE = 15 ! unit # for writing cmfd balance file + integer, parameter :: UNIT_PARTICLE = 16 ! unit # for writing particle restart + integer, parameter :: UNIT_OUTPUT = 17 ! unit # for writing output !============================================================================= ! CMFD CONSTANTS diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index d102648a8e..fc44489d00 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -17,18 +17,17 @@ module eigenvalue use random_lcg, only: prn, set_particle_seed, prn_skip use search, only: binary_search use source, only: get_source_particle - use state_point, only: write_state_point, replay_batch_history + use state_point, only: write_state_point use string, only: to_str use tally, only: synchronize_tallies, setup_active_usertallies, & reset_result -#ifdef HDF5 - use hdf5_interface, only: hdf5_write_state_point -#endif - private public :: run_eigenvalue + real(8) :: keff_generation ! Single-generation k on each processor + real(8) :: k_sum(2) = ZERO ! used to reduce sum and sum_sq + contains !=============================================================================== @@ -85,16 +84,7 @@ contains ! Accumulate time for transport call time_transport % stop() - ! Distribute fission bank across processors evenly - call time_bank % start() - call synchronize_bank() - call time_bank % stop() - - ! Calculate shannon entropy - if (entropy_on) call shannon_entropy() - - ! Write generation output - if (master .and. current_gen /= gen_per_batch) call print_generation() + call finalize_generation() end do GENERATION_LOOP @@ -147,14 +137,43 @@ contains subroutine initialize_generation() + ! set overall generation number + overall_gen = gen_per_batch*(current_batch - 1) + current_gen + ! Reset number of fission bank sites n_bank = 0 ! Count source sites if using uniform fission source weighting if (ufs) call count_source_for_ufs() + ! Store current value of tracklength k + keff_generation = global_tallies(K_TRACKLENGTH) % value + end subroutine initialize_generation +!=============================================================================== +! FINALIZE_GENERATION +!=============================================================================== + + subroutine finalize_generation() + + ! Distribute fission bank across processors evenly + call time_bank % start() + call synchronize_bank() + call time_bank % stop() + + ! Calculate shannon entropy + if (entropy_on) call shannon_entropy() + + ! Collect results and statistics + call calculate_generation_keff() + call calculate_average_keff() + + ! Write generation output + if (master .and. current_gen /= gen_per_batch) call print_generation() + + end subroutine finalize_generation + !=============================================================================== ! FINALIZE_BATCH handles synchronization and accumulation of tallies, ! calculation of Shannon entropy, getting single-batch estimate of keff, and @@ -168,8 +187,11 @@ contains call synchronize_tallies() call time_tallies % stop() - ! Collect results and statistics - call calculate_keff() + ! Reset global tally results + if (.not. active_batches) then + call reset_result(global_tallies) + n_realizations = 0 + end if ! Perform CMFD calculation if on if (cmfd_on) call execute_cmfd() @@ -183,11 +205,7 @@ contains if (master) call calculate_combined_keff() ! Create state point file -#ifdef HDF5 - call hdf5_write_state_point() -#else call write_state_point() -#endif end if if (master .and. current_batch == n_batches) then @@ -542,83 +560,59 @@ contains end subroutine shannon_entropy !=============================================================================== -! CALCULATE_KEFF calculates the single batch estimate of keff as well as the -! mean and standard deviation of the mean for active batches +! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto +! the master processor and normalizes them. This should work whether or not the +! no-reduce method is being used. !=============================================================================== - subroutine calculate_keff() + subroutine calculate_generation_keff() - real(8) :: temp(2) ! used to reduce sum and sum_sq - real(8) :: alpha ! significance level for CI - real(8) :: t_value ! t-value for confidence intervals - - message = "Calculate batch keff..." - call write_message(8) - - ! ========================================================================= - ! SINGLE-BATCH ESTIMATE OF K-EFFECTIVE + ! Get keff for this generation by subtracting off the starting value + keff_generation = global_tallies(K_TRACKLENGTH) % value - keff_generation #ifdef MPI - if (.not. reduce_tallies) then - ! Reduce value of k_batch if running in parallel - if (master) then - call MPI_REDUCE(MPI_IN_PLACE, k_batch(current_batch), 1, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) - else - ! Receive buffer not significant at other processors - call MPI_REDUCE(k_batch(current_batch), temp, 1, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) - end if - end if + ! Combine values across all processors + call MPI_ALLREDUCE(keff_generation, k_generation(overall_gen), 1, & + MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, mpi_err) +#else + k_generation(overall_gen) = keff_generation #endif ! Normalize single batch estimate of k - if (master) then - k_batch(current_batch) = k_batch(current_batch) / & - (n_particles * gen_per_batch) - end if + ! TODO: This should be normalized by total_weight, not by n_particles + k_generation(overall_gen) = k_generation(overall_gen) / n_particles - ! ========================================================================= - ! ACCUMULATED K-EFFECTIVE AND ITS VARIANCE + end subroutine calculate_generation_keff - if (reduce_tallies) then - ! In this case, global_tallies has already been reduced, so we don't - ! need to perform any more reductions and just take the values from - ! global_tallies directly +!=============================================================================== +! CALCULATE_AVERAGE_KEFF calculates the mean and standard deviation of the mean +! of k-effective during active generations and broadcasts the mean to all +! processors +!=============================================================================== - ! Sample mean of keff - keff = global_tallies(K_TRACKLENGTH) % sum / n_realizations + subroutine calculate_average_keff() - if (n_realizations > 1) then - if (confidence_intervals) then - ! Calculate t-value for confidence intervals - alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) - else - t_value = ONE - end if + integer :: n ! number of active generations + real(8) :: alpha ! significance level for CI + real(8) :: t_value ! t-value for confidence intervals + + ! Determine number of active generations + n = overall_gen - n_inactive*gen_per_batch + + if (n <= 0) then + ! For inactive generations, use current generation k as estimate for next + ! generation + keff = k_generation(overall_gen) - ! Standard deviation of the sample mean of k - keff_std = t_value * sqrt((global_tallies(K_TRACKLENGTH) % sum_sq / & - n_realizations - keff * keff) / (n_realizations - 1)) - end if else - ! In this case, no reduce was ever done on global_tallies. Thus, we need - ! to reduce the values in sum and sum^2 to get the sample mean and its - ! standard deviation + ! Sample mean of keff + k_sum(1) = k_sum(1) + k_generation(overall_gen) + k_sum(2) = k_sum(2) + k_generation(overall_gen)**2 -#ifdef MPI - call MPI_REDUCE(global_tallies(K_TRACKLENGTH) % sum, temp, 2, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) -#else - temp(1) = global_tallies(K_TRACKLENGTH) % sum - temp(2) = global_tallies(K_TRACKLENGTH) % sum_sq -#endif + ! Determine mean + keff = k_sum(1) / n - ! Sample mean of k - keff = temp(1) / n_realizations - - if (n_realizations > 1) then + if (n > 1) then if (confidence_intervals) then ! Calculate t-value for confidence intervals alpha = ONE - CONFIDENCE_LEVEL @@ -628,23 +622,11 @@ contains end if ! Standard deviation of the sample mean of k - keff_std = t_value * sqrt((temp(2)/n_realizations - keff*keff) / & - (n_realizations - 1)) + keff_std = t_value * sqrt((k_sum(2)/n - keff**2) / (n - 1)) end if end if -#ifdef MPI - ! Broadcast new keff value to all processors - call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) -#endif - - ! Reset global tally results - if (.not. active_batches) then - call reset_result(global_tallies) - n_realizations = 0 - end if - - end subroutine calculate_keff + end subroutine calculate_average_keff !=============================================================================== ! CALCULATE_COMBINED_KEFF calculates a minimum variance estimate of k-effective @@ -793,4 +775,37 @@ contains end subroutine count_source_for_ufs +!=============================================================================== +! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a +! batch using data read from a state point file +!=============================================================================== + + subroutine replay_batch_history + + ! Write message at beginning + if (current_batch == 1) then + message = "Replaying history from state point..." + call write_message(1) + end if + + do current_gen = 1, gen_per_batch + overall_gen = overall_gen + 1 + call calculate_average_keff() + + ! print out batch keff + if (current_gen < gen_per_batch) then + if (master) call print_generation() + else + if (master) call print_batch_keff() + end if + end do + + ! Write message at end + if (current_batch == restart_batch) then + message = "Resuming simulation..." + call write_message(1) + end if + + end subroutine replay_batch_history + end module eigenvalue diff --git a/src/finalize.F90 b/src/finalize.F90 index d601c2e592..bf0c752ffb 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -12,7 +12,7 @@ module finalize #endif #ifdef HDF5 - use hdf5_interface, only: hdf5_finalize + use hdf5_interface, only: h5tclose_f, h5close_f, hdf5_err #endif implicit none @@ -56,8 +56,12 @@ contains call free_memory() #ifdef HDF5 - ! Close HDF5 interface and release memory - call hdf5_finalize() + ! Release compound datatypes + call h5tclose_f(hdf5_tallyresult_t, hdf5_err) + call h5tclose_f(hdf5_bank_t, hdf5_err) + + ! Close FORTRAN interface. + call h5close_f(hdf5_err) #endif #ifdef MPI diff --git a/src/fixed_source.F90 b/src/fixed_source.F90 index 1d2b3bf00d..73c5291eb5 100644 --- a/src/fixed_source.F90 +++ b/src/fixed_source.F90 @@ -11,10 +11,6 @@ module fixed_source use string, only: to_str use tally, only: synchronize_tallies, setup_active_usertallies -#ifdef HDF5 - use hdf5_interface, only: hdf5_write_state_point -#endif - type(Bank), pointer :: source_site => null() contains @@ -115,11 +111,7 @@ contains ! Write out state point if it's been specified for this batch if (statepoint_batch % contains(current_batch)) then -#ifdef HDF5 - call hdf5_write_state_point() -#else call write_state_point() -#endif end if end subroutine finalize_batch diff --git a/src/global.F90 b/src/global.F90 index e0f0bd5bbf..5597ae2dfb 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -21,7 +21,7 @@ module global #endif #ifdef HDF5 - use hdf5 + use hdf5_interface, only: HID_T #endif implicit none @@ -156,7 +156,8 @@ module global integer :: n_active ! # of active batches integer :: gen_per_batch = 1 ! # of generations per batch integer :: current_batch = 0 ! current batch - integer :: current_gen = 0 ! current generation + integer :: current_gen = 0 ! current generation within a batch + integer :: overall_gen = 0 ! overall generation in the run ! External source type(ExtSource), target :: external_source @@ -172,7 +173,7 @@ module global integer(8) :: current_work ! index in source bank of current history simulated ! Temporary k-effective values - real(8), allocatable :: k_batch(:) ! batch estimates of k + real(8), allocatable :: k_generation(:) ! single-generation estimates of k real(8) :: keff = ONE ! average k over active batches real(8) :: keff_std ! standard deviation of average k real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption @@ -182,7 +183,7 @@ module global ! Shannon entropy logical :: entropy_on = .false. - real(8), allocatable :: entropy(:) ! shannon entropy at each batch + real(8), allocatable :: entropy(:) ! shannon entropy at each generation real(8), allocatable :: entropy_p(:,:,:,:) ! % of source sites in each cell type(StructuredMesh), pointer :: entropy_mesh @@ -244,7 +245,6 @@ module global integer(HID_T) :: hdf5_tallyresult_t ! Compound type for TallyResult integer(HID_T) :: hdf5_bank_t ! Compound type for Bank integer(HID_T) :: hdf5_integer8_t ! type for integer(8) - integer :: hdf5_err ! error flag #endif ! ============================================================================ diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 2aa5653fe8..471a4397cb 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -1,1306 +1,296 @@ module hdf5_interface - use ace_header, only: Reaction, UrrData - use bank_header, only: Bank - use constants - use endf, only: reaction_name - use error, only: fatal_error - use geometry_header, only: Cell, Surface, Universe, Lattice - use global - use material_header, only: Material - use output, only: write_message, time_stamp - use string, only: to_str - use tally_header, only: TallyObject - #ifdef HDF5 + use hdf5 use h5lt use, intrinsic :: ISO_C_BINDING + +#ifdef MPI + use mpi #endif implicit none -#ifdef HDF5 + integer :: hdf5_err ! HDF5 error code + integer :: hdf5_rank ! rank of data + integer(HID_T) :: dset ! data set handle + integer(HID_T) :: dspace ! data or file space handle + integer(HID_T) :: hdf5_fh ! HDF5 file handle + integer(HID_T) :: temp_group ! temporary HDF5 group handle + integer(HID_T) :: memspace ! data space handle for individual procs + integer(HID_T) :: plist ! property list handle + integer(HSIZE_T) :: dims1(1) ! dims type for 1-D array + integer(HSIZE_T) :: dims2(2) ! dims type for 2-D array + integer(HSIZE_T) :: dims3(3) ! dims type for 3-D array + type(c_ptr) :: f_ptr ! pointer to data + + ! Generic HDF5 write procedure interface + interface hdf5_write_data + module procedure hdf5_write_double + module procedure hdf5_write_double_1Darray + module procedure hdf5_write_integer + module procedure hdf5_write_integer_1Darray + module procedure hdf5_write_long + module procedure hdf5_write_string + end interface hdf5_write_data + + ! Generic HDF5 read procedure interface + interface hdf5_read_data + module procedure hdf5_read_double + module procedure hdf5_read_double_1Darray + module procedure hdf5_read_integer + module procedure hdf5_read_integer_1Darray + module procedure hdf5_read_long + module procedure hdf5_read_string + end interface hdf5_read_data contains !=============================================================================== -! HDF5_INITIALIZE +! HDF5_FILE_CREATE creates HDF5 file !=============================================================================== - subroutine hdf5_initialize() + subroutine hdf5_file_create(filename, file_id) - type(TallyResult), target :: tmp(2) ! temporary TallyResult - type(Bank), target :: tmpb(2) ! temporary Bank - integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals - integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates + character(*), intent(in) :: filename ! name of file + integer(HID_T), intent(inout) :: file_id ! file handle - ! Initialize FORTRAN interface. - call h5open_f(hdf5_err) + ! Create the file + call h5fcreate_f(trim(filename), H5F_ACC_TRUNC_F, file_id, hdf5_err) - ! Create the compound datatype for TallyResult - call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err) - call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err) - call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, 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) - - ! Determine type for integer(8) - hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND) - - end subroutine hdf5_initialize + end subroutine hdf5_file_create !=============================================================================== -! HDF5_FINALIZE +! HDF5_FILE_OPEN opens HDF5 file !=============================================================================== - subroutine hdf5_finalize() + subroutine hdf5_file_open(filename, file_id, mode) - ! Release compound datatypes - call h5tclose_f(hdf5_tallyresult_t, hdf5_err) + character(*), intent(in) :: filename ! name of file + character(*), intent(in) :: mode ! access mode to file + integer(HID_T), intent(inout) :: file_id ! file handle - ! Close FORTRAN interface. - call h5close_f(hdf5_err) + integer :: open_mode ! HDF5 open mode - end subroutine hdf5_finalize - -!=============================================================================== -! HDF5_WRITE_SUMMARY -!=============================================================================== - - subroutine hdf5_write_summary() - - character(MAX_FILE_LEN) :: filename = "summary.h5" - - ! Create a new file using default properties. - call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_output_file, hdf5_err) - - ! Write header information - call hdf5_write_header() - - ! Write eigenvalue information - if (run_mode == MODE_EIGENVALUE) then - ! Need to write integer(8)'s using double instead since there is no H5LT - ! call for making a dataset of type long - call hdf5_write_double(hdf5_output_file, "n_particles", real(n_particles,8)) - - ! Use H5LT interface to write n_batches, n_inactive, and n_active - call hdf5_write_integer(hdf5_output_file, "n_batches", n_batches) - call hdf5_write_integer(hdf5_output_file, "n_inactive", n_inactive) - call hdf5_write_integer(hdf5_output_file, "n_active", n_active) - call hdf5_write_integer(hdf5_output_file, "gen_per_batch", gen_per_batch) - - ! Add description of each variable - call h5ltset_attribute_string_f(hdf5_output_file, "n_particles", & - "description", "Number of particles per generation", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "n_batches", & - "description", "Total number of batches", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "n_inactive", & - "description", "Number of inactive batches", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "n_active", & - "description", "Number of active batches", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "gen_per_batch", & - "description", "Number of generations per batch", hdf5_err) + ! Determine access type + open_mode = H5F_ACC_RDONLY_F + if (trim(mode) == 'w') then + open_mode = H5F_ACC_RDWR_F end if - call hdf5_write_geometry() - call hdf5_write_materials() - call hdf5_write_nuclides() - if (n_tallies > 0) then - call hdf5_write_tallies() + ! Open file + call h5fopen_f(trim(filename), open_mode, file_id, hdf5_err) + + end subroutine hdf5_file_open + +!=============================================================================== +! HDF5_FILE_CLOSE closes HDF5 file +!=============================================================================== + + subroutine hdf5_file_close(file_id) + + integer(HID_T), intent(inout) :: file_id ! file handle + + ! Close the file + call h5fclose_f(file_id, hdf5_err) + + end subroutine hdf5_file_close + +#ifdef MPI + +!=============================================================================== +! HDF5_PARALLEL_FILE_CREATE creates HDF5 file with parallel I/O +!=============================================================================== + + subroutine hdf5_parallel_file_create(filename, file_id) + + character(*), intent(in) :: filename ! name of file + integer(HID_T), intent(inout) :: file_id ! file handle + + ! Setup file access property list with parallel I/O access + call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) + call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdf5_err) + + ! 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) + + end subroutine hdf5_parallel_file_create + +!=============================================================================== +! HDF5_PARALLEL_FILE_OPEN opens HDF5 file with parallel I/O +!=============================================================================== + + subroutine hdf5_parallel_file_open(filename, file_id, mode) + + character(*), intent(in) :: filename ! name of file + character(*), intent(in) :: mode ! access mode + integer(HID_T), intent(inout) :: file_id ! file handle + + integer :: open_mode ! HDF5 access mode + + ! Setup file access property list with parallel I/O access + call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) + call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdf5_err) + + ! Determine access type + open_mode = H5F_ACC_RDONLY_F + if (trim(mode) == 'w') then + open_mode = H5F_ACC_RDWR_F end if - ! Terminate access to the file. - call h5fclose_f(hdf5_output_file, hdf5_err) + ! Create the file collectively + call h5fopen_f(trim(filename), open_mode, file_id, hdf5_err, & + access_prp = plist) - end subroutine hdf5_write_summary + ! Close the property list + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_file_open + +#endif !=============================================================================== -! HDF5_WRITE_HEADER +! HDF5_OPEN_GROUP creates/opens HDF5 group to temp_group !=============================================================================== - subroutine hdf5_write_header() + subroutine hdf5_open_group(group) - ! Write version information - call hdf5_write_integer(hdf5_output_file, "version_major", VERSION_MAJOR) - call hdf5_write_integer(hdf5_output_file, "version_minor", VERSION_MINOR) - call hdf5_write_integer(hdf5_output_file, "version_release", VERSION_RELEASE) + character(*), intent(in) :: group ! name of group - ! Write current date and time - call h5ltmake_dataset_string_f(hdf5_output_file, "date_and_time", & - time_stamp(), hdf5_err) + logical :: status ! does the group exist - ! Write MPI information - call hdf5_write_integer(hdf5_output_file, "n_procs", n_procs) - call h5ltset_attribute_string_f(hdf5_output_file, "n_procs", & - "description", "Number of MPI processes", hdf5_err) + ! Check if group exists + call h5ltpath_valid_f(hdf5_fh, trim(group), .true., status, hdf5_err) - end subroutine hdf5_write_header - -!=============================================================================== -! HDF5_WRITE_GEOMETRY -!=============================================================================== - - subroutine hdf5_write_geometry() - - integer :: i, j, k, m - integer :: n_x, n_y, n_z - integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: dims3(3) - integer(HID_T) :: geometry_group - integer(HID_T) :: cell_group - integer(HID_T) :: surface_group - integer(HID_T) :: universe_group - integer(HID_T) :: lattice_group - integer(HID_T) :: temp_group - integer, allocatable :: lattice_universes(:,:,:) - type(Cell), pointer :: c => null() - type(Surface), pointer :: s => null() - type(Universe), pointer :: u => null() - type(Lattice), pointer :: lat => null() - - ! Create group for geometry - call h5gcreate_f(hdf5_output_file, "/geometry", geometry_group, hdf5_err) - - ! Use H5LT interface to write number of geometry objects - call hdf5_write_integer(geometry_group, "n_cells", n_cells) - call hdf5_write_integer(geometry_group, "n_surfaces", n_surfaces) - call hdf5_write_integer(geometry_group, "n_universes", n_universes) - call hdf5_write_integer(geometry_group, "n_lattices", n_lattices) - - ! ========================================================================== - ! WRITE INFORMATION ON CELLS - - call h5gcreate_f(geometry_group, "cells", cell_group, hdf5_err) - - ! Write information on each cell - do i = 1, n_cells - c => cells(i) - - ! Create group for i-th cell - call h5gcreate_f(cell_group, "cell " // trim(to_str(c % id)), & - temp_group, hdf5_err) - - ! Write universe for this cell - call hdf5_write_integer(temp_group, "universe", & - universes(c % universe) % id) - - ! Write information on what fills this cell - select case (c % type) - case (CELL_NORMAL) - call h5ltmake_dataset_string_f(temp_group, "fill_type", "normal", & - hdf5_err) - if (c % material == MATERIAL_VOID) then - call hdf5_write_integer(temp_group, "material", -1) - else - call hdf5_write_integer(temp_group, "material", & - materials(c % material) % id) - end if - case (CELL_FILL) - call h5ltmake_dataset_string_f(temp_group, "fill_type", "universe", & - hdf5_err) - call hdf5_write_integer(temp_group, "material", & - universes(c % fill) % id) - case (CELL_LATTICE) - call h5ltmake_dataset_string_f(temp_group, "fill_type", "lattice", & - hdf5_err) - call hdf5_write_integer(temp_group, "lattice", & - lattices(c % fill) % id) - end select - - ! Write list of bounding surfaces - if (c % n_surfaces > 0) then - dims(1) = c % n_surfaces - call h5ltmake_dataset_int_f(temp_group, "surfaces", 1, & - dims, c % surfaces, hdf5_err) - end if - - ! Close group for i-th cell - call h5gclose_f(temp_group, hdf5_err) - end do - - call h5gclose_f(cell_group, hdf5_err) - - ! ========================================================================== - ! WRITE INFORMATION ON SURFACES - - call h5gcreate_f(geometry_group, "surfaces", surface_group, hdf5_err) - - ! Write information on each surface - do i = 1, n_surfaces - s => surfaces(i) - - ! Create group for i-th surface - call h5gcreate_f(surface_group, "surface " // trim(to_str(s % id)), & - temp_group, hdf5_err) - - ! Write surface type - select case (s % type) - case (SURF_PX) - call h5ltmake_dataset_string_f(temp_group, "type", "X Plane", hdf5_err) - case (SURF_PY) - call h5ltmake_dataset_string_f(temp_group, "type", "Y Plane", hdf5_err) - case (SURF_PZ) - call h5ltmake_dataset_string_f(temp_group, "type", "Z Plane", hdf5_err) - case (SURF_PLANE) - call h5ltmake_dataset_string_f(temp_group, "type", "Plane", hdf5_err) - case (SURF_CYL_X) - call h5ltmake_dataset_string_f(temp_group, "type", "X Cylinder", hdf5_err) - case (SURF_CYL_Y) - call h5ltmake_dataset_string_f(temp_group, "type", "Y Cylinder", hdf5_err) - case (SURF_CYL_Z) - call h5ltmake_dataset_string_f(temp_group, "type", "Z Cylinder", hdf5_err) - case (SURF_SPHERE) - call h5ltmake_dataset_string_f(temp_group, "type", "Sphere", hdf5_err) - case (SURF_CONE_X) - call h5ltmake_dataset_string_f(temp_group, "type", "X Cone", hdf5_err) - case (SURF_CONE_Y) - call h5ltmake_dataset_string_f(temp_group, "type", "Y Cone", hdf5_err) - case (SURF_CONE_Z) - call h5ltmake_dataset_string_f(temp_group, "type", "Z Cone", hdf5_err) - end select - - ! Write coefficients for surface - dims(1) = size(s % coeffs) - call h5ltmake_dataset_double_f(temp_group, "coefficients", 1, dims, & - s % coeffs, hdf5_err) - - ! Write positive neighbors - if (allocated(s % neighbor_pos)) then - dims(1) = size(s % neighbor_pos) - call h5ltmake_dataset_int_f(temp_group, "neighbors_positive", 1, dims, & - s % neighbor_pos, hdf5_err) - end if - - ! Write negative neighbors - if (allocated(s % neighbor_neg)) then - dims(1) = size(s % neighbor_neg) - call h5ltmake_dataset_int_f(temp_group, "neighbors_negative", 1, dims, & - s % neighbor_neg, hdf5_err) - end if - - ! Write boundary condition - select case (s % bc) - case (BC_TRANSMIT) - call h5ltmake_dataset_string_f(temp_group, "boundary_condition", & - "transmission", hdf5_err) - case (BC_VACUUM) - call h5ltmake_dataset_string_f(temp_group, "boundary_condition", & - "vacuum", hdf5_err) - case (BC_REFLECT) - call h5ltmake_dataset_string_f(temp_group, "boundary_condition", & - "reflective", hdf5_err) - case (BC_PERIODIC) - call h5ltmake_dataset_string_f(temp_group, "boundary_condition", & - "periodic", hdf5_err) - end select - - ! Close group for i-th surface - call h5gclose_f(temp_group, hdf5_err) - end do - - call h5gclose_f(surface_group, hdf5_err) - - ! ========================================================================== - ! WRITE INFORMATION ON UNIVERSES - - call h5gcreate_f(geometry_group, "universes", universe_group, hdf5_err) - - ! Write information on each universe - do i = 1, n_universes - u => universes(i) - - ! Create group for i-th universe - call h5gcreate_f(universe_group, "universe " // trim(to_str(u % id)), & - temp_group, hdf5_err) - - ! Write list of cells in this universe - if (u % n_cells > 0) then - dims(1) = u % n_cells - call h5ltmake_dataset_int_f(temp_group, "cells", 1, dims, & - u % cells, hdf5_err) - end if - - ! Close group for i-th universe - call h5gclose_f(temp_group, hdf5_err) - end do - - call h5gclose_f(universe_group, hdf5_err) - - ! ========================================================================== - ! WRITE INFORMATION ON LATTICES - - call h5gcreate_f(geometry_group, "lattices", lattice_group, hdf5_err) - - ! Write information on each lattice - do i = 1, n_lattices - lat => lattices(i) - - ! Create group for i-th lattice - call h5gcreate_f(lattice_group, "lattice " // trim(to_str(lat % id)), & - temp_group, hdf5_err) - - ! Write lattice type - select case(lat % type) - case (LATTICE_RECT) - call h5ltmake_dataset_string_f(temp_group, "type", "rectangular", hdf5_err) - case (LATTICE_HEX) - call h5ltmake_dataset_string_f(temp_group, "type", "hexagonal", hdf5_err) - end select - - ! Write lattice dimensions, lower left corner, and width of element - dims(1) = lat % n_dimension - call h5ltmake_dataset_int_f(temp_group, "dimension", 1, dims, & - lat % dimension, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "lower_left", 1, dims, & - lat % lower_left, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "width", 1, dims, & - lat % width, hdf5_err) - - ! Determine dimensions of lattice - n_x = lat % dimension(1) - n_y = lat % dimension(2) - if (lat % n_dimension == 3) then - n_z = lat % dimension(3) - else - n_z = 1 - end if - - ! Write lattice universes - allocate(lattice_universes(n_x, n_y, n_z)) - do j = 1, n_x - do k = 1, n_y - do m = 1, n_z - lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id - end do - end do - end do - dims3 = (/ n_x, n_y, n_z /) - call h5ltmake_dataset_int_f(temp_group, "universes", 3, dims3, & - lattice_universes, hdf5_err) - deallocate(lattice_universes) - - ! Close group for i-th lattice - call h5gclose_f(temp_group, hdf5_err) - end do - - call h5gclose_f(lattice_group, hdf5_err) - - ! Close geometry group - call h5gclose_f(geometry_group, hdf5_err) - - end subroutine hdf5_write_geometry - -!=============================================================================== -! HDF5_WRITE_MATERIALS -!=============================================================================== - - subroutine hdf5_write_materials() - - integer :: i - integer :: j - integer(HSIZE_T) :: dims(1) - integer(HID_T) :: materials_group - integer(HID_T) :: temp_group - integer, allocatable :: zaids(:) - type(Material), pointer :: m => null() - - ! Create group for materials - call h5gcreate_f(hdf5_output_file, "/materials", materials_group, hdf5_err) - - ! Use H5LT interface to write number of materials - call hdf5_write_integer(materials_group, "n_materials", n_materials) - - ! Write information on each material - do i = 1, n_materials - m => materials(i) - - ! Create group for i-th universe - call h5gcreate_f(materials_group, "material " // trim(to_str(m % id)), & - temp_group, hdf5_err) - - ! Write atom density with units - call hdf5_write_double(temp_group, "atom_density", m % density) - call h5ltset_attribute_string_f(temp_group, "atom_density", & - "units", "atom/barn-cm", hdf5_err) - - ! Copy ZAID for each nuclide to temporary array - allocate(zaids(m % n_nuclides)) - do j = 1, m % n_nuclides - zaids(j) = nuclides(m % nuclide(j)) % zaid - end do - - ! Write temporary array to 'nuclides' - dims(1) = m % n_nuclides - call h5ltmake_dataset_int_f(temp_group, "nuclides", 1, & - dims, zaids, hdf5_err) - - ! Deallocate temporary array - deallocate(zaids) - - ! Write atom densities - call h5ltmake_dataset_double_f(temp_group, "nuclide_densities", 1, & - dims, m % atom_density, hdf5_err) - - ! Write S(a,b) information if present - if (m % n_sab > 0) then - dims(1) = m % n_sab - call h5ltmake_dataset_int_f(temp_group, "i_sab_nuclides", 1, & - dims, m % i_sab_nuclides, hdf5_err) - call h5ltmake_dataset_int_f(temp_group, "i_sab_tables", 1, & - dims, m % i_sab_tables, hdf5_err) - end if - - ! Close group for i-th material - call h5gclose_f(temp_group, hdf5_err) - end do - - ! Close materials group - call h5gclose_f(materials_group, hdf5_err) - - end subroutine hdf5_write_materials - -!=============================================================================== -! HDF5_WRITE_TALLIES -!=============================================================================== - - subroutine hdf5_write_tallies() - - integer :: i, j - integer(HSIZE_T) :: dims(1) - integer(HID_T) :: tallies_group - integer(HID_T) :: temp_group - integer(HID_T) :: filter_group ! group for i-th filter - integer, allocatable :: temp_array(:) ! nuclide bin array - type(TallyObject), pointer :: t => null() - - ! Create group for tallies - call h5gcreate_f(hdf5_output_file, "tallies", tallies_group, hdf5_err) - - ! Write total number of meshes - call hdf5_write_integer(tallies_group, "n_meshes", n_meshes) - - ! Write information for meshes - MESH_LOOP: do i = 1, n_meshes - ! Create temporary group for each mesh - call h5gcreate_f(tallies_group, "mesh" // to_str(i), & - temp_group, hdf5_err) - - ! Write type and number of dimensions - call hdf5_write_integer(temp_group, "type", meshes(i) % type) - call hdf5_write_integer(temp_group, "n_dimension", & - meshes(i) % n_dimension) - - ! Write mesh information - dims(1) = meshes(i) % n_dimension - call h5ltmake_dataset_int_f(temp_group, "dimension", 1, & - dims, meshes(i) % dimension, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "lower_left", 1, & - dims, meshes(i) % lower_left, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "upper_right", 1, & - dims, meshes(i) % upper_right, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "width", 1, & - dims, meshes(i) % width, hdf5_err) - - ! Close temporary group for mesh - call h5gclose_f(temp_group, hdf5_err) - end do MESH_LOOP - - ! Write number of tallies - call hdf5_write_integer(tallies_group, "n_tallies", n_tallies) - - TALLY_METADATA: do i = 1, n_tallies - ! Get pointer to tally - t => tallies(i) - - ! Create group for this tally - call h5gcreate_f(tallies_group, "tally" // to_str(i), & - temp_group, hdf5_err) - - ! Write size of each tally - call hdf5_write_integer(temp_group, "total_score_bins", & - t % total_score_bins) - call hdf5_write_integer(temp_group, "total_filter_bins", & - t % total_filter_bins) - - ! Write number of filters - call hdf5_write_integer(temp_group, "n_filters", t % n_filters) - - FILTER_LOOP: do j = 1, t % n_filters - ! Create filter group - call h5gcreate_f(temp_group, "filter" // to_str(j), filter_group, & - hdf5_err) - - ! Write type of filter - call hdf5_write_integer(filter_group, "type", t % filters(j) % type) - - ! Write number of bins for this filter - call hdf5_write_integer(filter_group, "n_bins", t % filters(j) % n_bins) - - ! Write filter bins - if (t % filters(j) % type == FILTER_ENERGYIN .or. & - t % filters(j) % type == FILTER_ENERGYOUT) then - dims(1) = size(t % filters(j) % real_bins) - call h5ltmake_dataset_double_f(filter_group, "bins", 1, & - dims, t % filters(j) % real_bins, hdf5_err) - else - dims(1) = size(t % filters(j) % int_bins) - call h5ltmake_dataset_int_f(filter_group, "bins", 1, & - dims, t % filters(j) % int_bins, hdf5_err) - end if - - ! Write name of type - select case (t % filters(j) % type) - case(FILTER_UNIVERSE) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "universe", hdf5_err) - case(FILTER_MATERIAL) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "material", hdf5_err) - case(FILTER_CELL) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "cell", hdf5_err) - case(FILTER_CELLBORN) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "cellborn", hdf5_err) - case(FILTER_SURFACE) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "surface", hdf5_err) - case(FILTER_MESH) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "mesh", hdf5_err) - case(FILTER_ENERGYIN) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "energy", hdf5_err) - case(FILTER_ENERGYOUT) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "energyout", hdf5_err) - end select - - ! Close group for this filter - call h5gclose_f(filter_group, hdf5_err) - end do FILTER_LOOP - - ! Write number of nuclide bins - call hdf5_write_integer(temp_group, "n_nuclide_bins", & - t % n_nuclide_bins) - - - ! Create temporary array for nuclide bins - allocate(temp_array(t % n_nuclide_bins)) - NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins - if (t % nuclide_bins(j) > 0) then - temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid - else - temp_array(j) = t % nuclide_bins(j) - end if - end do NUCLIDE_LOOP - - ! Write and deallocate nuclide bins - dims(1) = t % n_nuclide_bins - call h5ltmake_dataset_int_f(temp_group, "nuclide_bins", 1, & - dims, temp_array, hdf5_err) - deallocate(temp_array) - - ! Write number of score bins - call hdf5_write_integer(temp_group, "n_score_bins", & - t % n_score_bins) - dims(1) = t % n_score_bins - call h5ltmake_dataset_int_f(temp_group, "score_bins", 1, & - dims, t % score_bins, hdf5_err) - - ! Close tally group - call h5gclose_f(temp_group, hdf5_err) - end do TALLY_METADATA - - ! Close tallies group - call h5gclose_f(tallies_group, hdf5_err) - - end subroutine hdf5_write_tallies - -!=============================================================================== -! HDF5_WRITE_NUCLIDES -!=============================================================================== - - subroutine hdf5_write_nuclides() - - integer :: i, j - integer :: size_total - integer :: size_xs - integer :: size_angle - integer :: size_energy - integer(HID_T) :: group - integer(HID_T) :: nuclide_group - integer(HID_T) :: reactions_group - integer(HID_T) :: rxn_group - type(Nuclide), pointer :: nuc => null() - type(Reaction), pointer :: rxn => null() - type(UrrData), pointer :: urr => null() - - ! Create group for nuclides - call h5gcreate_f(hdf5_output_file, "/nuclides", group, hdf5_err) - - ! Use H5LT interface to write number of nuclides - call hdf5_write_integer(group, "n_nuclides", n_nuclides_total) - - ! Write information on each nuclide - do i = 1, n_nuclides_total - nuc => nuclides(i) - - ! Determine size of cross-sections - size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8 - size_total = size_xs - - ! Create group for i-th nuclide - call h5gcreate_f(group, trim(nuc % name), nuclide_group, hdf5_err) - - ! Write some basic attributes - call hdf5_write_integer(nuclide_group, "zaid", nuc % zaid) - call hdf5_write_double(nuclide_group, "awr", nuc % awr) - call hdf5_write_double(nuclide_group, "kT", nuc % kT) - call hdf5_write_integer(nuclide_group, "n_grid", nuc % n_grid) - call hdf5_write_integer(nuclide_group, "n_reactions", nuc % n_reaction) - call hdf5_write_integer(nuclide_group, "n_fission", nuc % n_fission) - call hdf5_write_integer(nuclide_group, "size_xs", size_xs) - - ! ======================================================================= - ! WRITE INFORMATION ON EACH REACTION - - ! Create overall group - call h5gcreate_f(nuclide_group, "reactions", reactions_group, hdf5_err) - - do j = 1, nuc % n_reaction - ! Information on each reaction - rxn => nuc % reactions(j) - - ! Determine size of angle distribution - if (rxn % has_angle_dist) then - size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 - else - size_angle = 0 - end if - - ! Determine size of energy distribution - if (rxn % has_energy_dist) then - size_energy = size(rxn % edist % data) * 8 - else - size_energy = 0 - end if - - ! Create reaction group - call h5gcreate_f(reactions_group, reaction_name(rxn % MT), & - rxn_group, hdf5_err) - - ! Write information on reaction - call hdf5_write_double(rxn_group, "Q_value", rxn % Q_value) - call hdf5_write_integer(rxn_group, "multiplicity", rxn % multiplicity) - call hdf5_write_integer(rxn_group, "threshold", rxn % threshold) - call hdf5_write_integer(rxn_group, "size_angle", size_angle) - call hdf5_write_integer(rxn_group, "size_energy", size_energy) - - call h5gclose_f(rxn_group, hdf5_err) - - ! Accumulate data size - size_total = size_total + size_angle + size_energy - end do - - ! Close overall group for reactions - call h5gclose_f(reactions_group, hdf5_err) - - ! ======================================================================= - ! WRITE INFORMATION ON URR PROBABILITY TABLES - - if (nuc % urr_present) then - urr => nuc % urr_data - call hdf5_write_integer(nuclide_group, "urr_n_energy", urr % n_energy) - call hdf5_write_integer(nuclide_group, "urr_n_prob", urr % n_prob) - call hdf5_write_integer(nuclide_group, "urr_interp", urr % interp) - call hdf5_write_integer(nuclide_group, "urr_inelastic", urr % inelastic_flag) - call hdf5_write_integer(nuclide_group, "urr_absorption", urr % absorption_flag) - call hdf5_write_double(nuclide_group, "urr_min_E", urr % energy(1)) - call hdf5_write_double(nuclide_group, "urr_max_E", urr % energy(urr % n_energy)) - end if - - ! Write total memory used - call hdf5_write_integer(nuclide_group, "size_total", size_total) - - ! Close group for i-th nuclide - call h5gclose_f(nuclide_group, hdf5_err) - end do - - ! Close group for nuclides - call h5gclose_f(group, hdf5_err) - - end subroutine hdf5_write_nuclides - -!=============================================================================== -! HDF5_WRITE_TIMING -!=============================================================================== - - subroutine hdf5_write_timing() - - integer(HID_T) :: timing_group - integer(8) :: total_particles - real(8) :: speed - - ! Create group for timing - call h5gcreate_f(hdf5_output_file, "/timing", timing_group, hdf5_err) - - ! Write timing data - call hdf5_write_double(timing_group, "time_initialize", time_initialize % elapsed) - call hdf5_write_double(timing_group, "time_read_xs", time_read_xs % elapsed) - call hdf5_write_double(timing_group, "time_unionize", time_unionize % elapsed) - call hdf5_write_double(timing_group, "time_transport", time_transport % elapsed) - call hdf5_write_double(timing_group, "time_bank", time_bank % elapsed) - call hdf5_write_double(timing_group, "time_bank_sample", time_bank_sample % elapsed) - call hdf5_write_double(timing_group, "time_bank_sendrecv", time_bank_sendrecv % elapsed) - call hdf5_write_double(timing_group, "time_tallies", time_tallies % elapsed) - call hdf5_write_double(timing_group, "time_inactive", time_inactive % elapsed) - call hdf5_write_double(timing_group, "time_active", time_active % elapsed) - call hdf5_write_double(timing_group, "time_finalize", time_finalize % elapsed) - call hdf5_write_double(timing_group, "time_total", time_total % elapsed) - - ! Add descriptions to timing data - call h5ltset_attribute_string_f(timing_group, "time_initialize", & - "description", "Total time elapsed for initialization (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_read_xs", & - "description", "Time reading cross-section libraries (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_unionize", & - "description", "Time unionizing energy grid (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_transport", & - "description", "Time in transport only (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_bank", & - "description", "Total time synchronizing fission bank (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_bank_sample", & - "description", "Time between generations sampling source sites (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_bank_sendrecv", & - "description", "Time between generations SEND/RECVing source sites (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_tallies", & - "description", "Time between batches accumulating tallies (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_inactive", & - "description", "Total time in inactive batches (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_active", & - "description", "Total time in active batches (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_finalize", & - "description", "Total time for finalization (s)", hdf5_err) - call h5ltset_attribute_string_f(timing_group, "time_total", & - "description", "Total time elapsed (s)", hdf5_err) - - ! Write calculation rate - total_particles = n_particles * n_batches * gen_per_batch - speed = real(total_particles) / (time_inactive % elapsed + & - time_active % elapsed) - call hdf5_write_double(timing_group, "neutrons_per_second", speed) - - ! Close timing group - call h5gclose_f(timing_group, hdf5_err) - - end subroutine hdf5_write_timing - -!=============================================================================== -! HDF5_WRITE_STATE_POINT -!=============================================================================== - - subroutine hdf5_write_state_point() - - integer :: i, j ! loop indices - integer(HSIZE_T) :: dims(1) ! dimensions of 1D arrays - integer(HSIZE_T) :: dims2(2) ! dimensions of 2D arrays - integer(HSIZE_T) :: dims4(4) ! dimensions of 4D arrays - integer(HID_T) :: hdf5_state_point ! identifier for state point file - integer(HID_T) :: tallies_group ! "tallies" group - integer(HID_T) :: temp_group ! group for i-th tally or mesh - integer(HID_T) :: filter_group ! group for i-th filter - integer(HID_T) :: cmfd_group ! group for cmfd output - integer(HID_T) :: dspace ! identifier for dataspace - integer(HID_T) :: dset ! identifier for dataset - integer, allocatable :: temp_array(:) ! nuclide bin array - character(MAX_FILE_LEN) :: filename ! Path to state point - type(c_ptr) :: f_ptr ! Pointer for h5dwrite - type(TallyObject), pointer :: t => null() - - filename = trim(path_output) // 'statepoint.' // & - trim(to_str(current_batch)) // '.h5' - - ! Write message - message = "Creating HDF5 state point " // trim(filename) // "..." - call write_message(1) - - ! Only master process should continue from here - if (.not. master) return - - ! Create a new state point file using default properties. - call h5fcreate_f(filename, H5F_ACC_TRUNC_F, hdf5_state_point, & - hdf5_err) - - ! Write revision number for state point file - call hdf5_write_integer(hdf5_state_point, "revision_statepoint", & - REVISION_STATEPOINT) - - ! Write OpenMC version - call hdf5_write_integer(hdf5_state_point, "version_major", VERSION_MAJOR) - call hdf5_write_integer(hdf5_state_point, "version_minor", VERSION_MINOR) - call hdf5_write_integer(hdf5_state_point, "version_release", VERSION_RELEASE) - - ! Write current date and time - call h5ltmake_dataset_string_f(hdf5_state_point, "date_and_time", & - time_stamp(), hdf5_err) - - ! Write path to input - call h5ltmake_dataset_string_f(hdf5_state_point, "path", & - path_input, hdf5_err) - - ! Write out random number seed - call hdf5_write_long(hdf5_state_point, "seed", seed) - - ! Write run information - call hdf5_write_integer(hdf5_state_point, "run_mode", run_mode) - call hdf5_write_long(hdf5_state_point, "n_particles", n_particles) - call hdf5_write_integer(hdf5_state_point, "n_batches", n_batches) - - ! Write out current batch number - call hdf5_write_integer(hdf5_state_point, "current_batch", current_batch) - - ! Write out information for eigenvalue run - if (run_mode == MODE_EIGENVALUE) then - call hdf5_write_integer(hdf5_state_point, "n_inactive", n_inactive) - call hdf5_write_integer(hdf5_state_point, "gen_per_batch", gen_per_batch) - - ! Write out keff and entropy - dims(1) = current_batch - call h5ltmake_dataset_double_f(hdf5_state_point, "k_batch", 1, & - dims, k_batch, hdf5_err) - dims(1) = current_batch*gen_per_batch - call h5ltmake_dataset_double_f(hdf5_state_point, "entropy", 1, & - dims, entropy, hdf5_err) - call hdf5_write_double(hdf5_state_point, "k_col_abs", k_col_abs) - call hdf5_write_double(hdf5_state_point, "k_col_tra", k_col_tra) - call hdf5_write_double(hdf5_state_point, "k_abs_tra", k_abs_tra) - dims(1) = 2 - call h5ltmake_dataset_double_f(hdf5_state_point, "k_combined", 1, & - dims, k_combined, hdf5_err) - end if - - ! Create group for tallies - call h5gcreate_f(hdf5_state_point, "tallies", tallies_group, hdf5_err) - - ! Write total number of meshes - call hdf5_write_integer(tallies_group, "n_meshes", n_meshes) - - ! Write information for meshes - MESH_LOOP: do i = 1, n_meshes - ! Create temporary group for each mesh - call h5gcreate_f(tallies_group, "mesh" // to_str(i), & - temp_group, hdf5_err) - - ! Write id, type, and number of dimensions - call hdf5_write_integer(temp_group, "id", meshes(i) % id) - call hdf5_write_integer(temp_group, "type", meshes(i) % type) - call hdf5_write_integer(temp_group, "n_dimension", & - meshes(i) % n_dimension) - - ! Write mesh information - dims(1) = meshes(i) % n_dimension - call h5ltmake_dataset_int_f(temp_group, "dimension", 1, & - dims, meshes(i) % dimension, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "lower_left", 1, & - dims, meshes(i) % lower_left, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "upper_right", 1, & - dims, meshes(i) % upper_right, hdf5_err) - call h5ltmake_dataset_double_f(temp_group, "width", 1, & - dims, meshes(i) % width, hdf5_err) - - ! Close temporary group for mesh - call h5gclose_f(temp_group, hdf5_err) - end do MESH_LOOP - - ! Write number of tallies - call hdf5_write_integer(tallies_group, "n_tallies", n_tallies) - - TALLY_METADATA: do i = 1, n_tallies - ! Get pointer to tally - t => tallies(i) - - ! Create group for this tally - call h5gcreate_f(tallies_group, "tally" // to_str(i), & - temp_group, hdf5_err) - - ! Write id - call hdf5_write_integer(temp_group, "id", t % id) - - ! Write number of realizations - call hdf5_write_integer(temp_group, "n_realizations", & - t % n_realizations) - - ! Write size of each tally - call hdf5_write_integer(temp_group, "total_score_bins", & - t % total_score_bins) - call hdf5_write_integer(temp_group, "total_filter_bins", & - t % total_filter_bins) - - ! Write number of filters - call hdf5_write_integer(temp_group, "n_filters", t % n_filters) - - FILTER_LOOP: do j = 1, t % n_filters - ! Create filter group - call h5gcreate_f(temp_group, "filter" // to_str(j), filter_group, & - hdf5_err) - - ! Write type of filter - call hdf5_write_integer(filter_group, "type", t % filters(j) % type) - - ! Write number of bins for this filter - call hdf5_write_integer(filter_group, "n_bins", t % filters(j) % n_bins) - - ! Write filter bins - if (t % filters(j) % type == FILTER_ENERGYIN .or. & - t % filters(j) % type == FILTER_ENERGYOUT) then - dims(1) = size(t % filters(j) % real_bins) - call h5ltmake_dataset_double_f(filter_group, "bins", 1, & - dims, t % filters(j) % real_bins, hdf5_err) - else - dims(1) = size(t % filters(j) % int_bins) - call h5ltmake_dataset_int_f(filter_group, "bins", 1, & - dims, t % filters(j) % int_bins, hdf5_err) - end if - - ! Write name of type - select case (t % filters(j) % type) - case(FILTER_UNIVERSE) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "universe", hdf5_err) - case(FILTER_MATERIAL) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "material", hdf5_err) - case(FILTER_CELL) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "cell", hdf5_err) - case(FILTER_CELLBORN) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "cellborn", hdf5_err) - case(FILTER_SURFACE) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "surface", hdf5_err) - case(FILTER_MESH) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "mesh", hdf5_err) - case(FILTER_ENERGYIN) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "energy", hdf5_err) - case(FILTER_ENERGYOUT) - call h5ltmake_dataset_string_f(filter_group, "type_name", & - "energyout", hdf5_err) - end select - - ! Close group for this filter - call h5gclose_f(filter_group, hdf5_err) - end do FILTER_LOOP - - ! Write number of nuclide bins - call hdf5_write_integer(temp_group, "n_nuclide_bins", & - t % n_nuclide_bins) - - - ! Create temporary array for nuclide bins - allocate(temp_array(t % n_nuclide_bins)) - NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins - if (t % nuclide_bins(j) > 0) then - temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid - else - temp_array(j) = t % nuclide_bins(j) - end if - end do NUCLIDE_LOOP - - ! Write and deallocate nuclide bins - dims(1) = t % n_nuclide_bins - call h5ltmake_dataset_int_f(temp_group, "nuclide_bins", 1, & - dims, temp_array, hdf5_err) - deallocate(temp_array) - - ! Write number of score bins - call hdf5_write_integer(temp_group, "n_score_bins", & - t % n_score_bins) - - ! Write score bins and scattering order - dims(1) = t % n_score_bins - call h5ltmake_dataset_int_f(temp_group, "score_bins", 1, & - dims, t % score_bins, hdf5_err) - call h5ltmake_dataset_int_f(temp_group, "scatt_order", 1, & - dims, t % scatt_order, hdf5_err) - - ! Write number of user score bins - call hdf5_write_integer(temp_group, "n_user_score_bins", & - t % n_user_score_bins) - - ! Close tally group - call h5gclose_f(temp_group, hdf5_err) - end do TALLY_METADATA - - ! Write number of realizations for global tallies - call hdf5_write_integer(hdf5_state_point, "n_realizations", n_realizations) - - ! Write out global tallies sum and sum_sq - call hdf5_write_integer(hdf5_state_point, "n_global_tallies", & - N_GLOBAL_TALLIES) - - ! Write global tallies - dims(1) = N_GLOBAL_TALLIES - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyresult_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(global_tallies(1)) - CALL h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - - if (tallies_on) then - ! Indicate that tallies are on - call hdf5_write_integer(tallies_group, "tallies_present", 1) - - ! Write tally sum and sum_sq - TALLY_RESULTS: do i = 1, n_tallies - ! Get pointer to tally - t => tallies(i) - - ! Open group for the i-th tally - call h5gopen_f(tallies_group, "tally" // to_str(i), & - temp_group, hdf5_err) - - ! Write sum and sum_sq for each bin - dims2 = shape(t % results) - call h5screate_simple_f(2, dims2, dspace, hdf5_err) - call h5dcreate_f(temp_group, "results", hdf5_tallyresult_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(t % results(1, 1)) - CALL h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - - ! Close group for the i-th tally - call h5gclose_f(temp_group, hdf5_err) - end do TALLY_RESULTS + ! Either create or open group + if (status) then + call h5gopen_f(hdf5_fh, trim(group), temp_group, hdf5_err) else - ! Indicate that tallies are off - call hdf5_write_integer(tallies_group, "tallies_present", 0) + call h5gcreate_f(hdf5_fh, trim(group), temp_group, hdf5_err) end if - ! Close tallies group - call h5gclose_f(tallies_group, hdf5_err) - - ! TODO: Use parallel HDF5 to write source bank - - ! Write source bank - if (source_write) then - dims(1) = work - call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(hdf5_state_point, "source_bank", hdf5_bank_t, & - dspace, dset, hdf5_err) - f_ptr = c_loc(source_bank(1)) - CALL h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - call h5dclose_f(dset, hdf5_err) - call h5sclose_f(dspace, hdf5_err) - end if - - ! Write out CMFD info if active - if (cmfd_on) then - - ! Create CMFD group - call h5gcreate_f(hdf5_state_point, "cmfd", cmfd_group, hdf5_err) - - ! write out openmc source - dims4 = shape(cmfd % openmc_src) - call h5ltmake_dataset_double_f(cmfd_group, "openmc_src", 4, & - dims4, cmfd % openmc_src, hdf5_err) - - ! write out cmfd source - dims4 = shape(cmfd % cmfd_src) - call h5ltmake_dataset_double_f(cmfd_group, "cmfd_src", 4, & - dims4, cmfd % cmfd_src, hdf5_err) - - ! write out keff - call hdf5_write_double(cmfd_group, "cmfd_keff", cmfd % keff) - - ! Close CMFD group - call h5gclose_f(tallies_group, hdf5_err) - - end if - - ! Close HDF5 state point file - call h5fclose_f(hdf5_state_point, hdf5_err) - - end subroutine hdf5_write_state_point + end subroutine hdf5_open_group !=============================================================================== -! HDF5_LOAD_STATE_POINT +! HDF5_CLOSE_GROUP closes HDF5 temp_group !=============================================================================== - subroutine hdf5_load_state_point() + subroutine hdf5_close_group() - integer :: i ! loop index - integer :: mode ! specified run mode - integer :: temp ! statepoint revision for comparison - integer(HID_T) :: hdf5_state_point ! identifier for state point file - integer(HID_T) :: tally_group ! identifier for tally group - integer(HID_T) :: dset ! identifier for dataset - integer(HSIZE_T) :: dims(1) ! dimensions of 1D arrays - type(c_ptr) :: f_ptr ! c pointer for h5dread + ! Close the group + call h5gclose_f(temp_group, hdf5_err) - ! Write message - message = "Loading HDF5 state point " // trim(path_state_point) // "..." - call write_message(1) - - ! Load state point file - call h5fopen_f(path_state_point, H5F_ACC_RDONLY_F, & - hdf5_state_point, hdf5_err) - - ! Raad revision number for state point file and make sure it matches with - ! current version - call hdf5_read_integer(hdf5_state_point, "revision_statepoint", temp) - if (temp /= REVISION_STATEPOINT) then - message = "State point binary version does not match current version " & - // "in OpenMC." - call fatal_error() - end if - - ! Read and overwrite random number seed - call hdf5_read_long(hdf5_state_point, "seed", seed) - - ! Read and overwrite run information - call hdf5_read_integer(hdf5_state_point, "run_mode", mode) - call hdf5_read_long(hdf5_state_point, "n_particles", n_particles) - call hdf5_read_integer(hdf5_state_point, "n_batches", n_batches) - - ! Read batch number to restart at - call hdf5_read_integer(hdf5_state_point, "current_batch", restart_batch) - - ! Read information specific to eigenvalue run - if (mode == MODE_EIGENVALUE) then - call hdf5_read_integer(hdf5_state_point, "n_inactive", n_inactive) - call hdf5_read_integer(hdf5_state_point, "gen_per_batch", gen_per_batch) - - dims(1) = restart_batch - call h5ltread_dataset_double_f(hdf5_state_point, "k_batch", & - k_batch(1:restart_batch), dims, hdf5_err) - dims(1) = restart_batch*gen_per_batch - call h5ltread_dataset_double_f(hdf5_state_point, "entropy", & - entropy(1:restart_batch*gen_per_batch), dims, hdf5_err) - call hdf5_read_double(hdf5_state_point, "k_col_abs", k_col_abs) - call hdf5_read_double(hdf5_state_point, "k_col_tra", k_col_tra) - call hdf5_read_double(hdf5_state_point, "k_abs_tra", k_abs_tra) - end if - - ! Read number of realizations for global tallies - call hdf5_read_integer(hdf5_state_point, "n_realizations", n_realizations) - - if (master) then - ! Open global tallies dataset - call h5dopen_f(hdf5_state_point, "global_tallies", dset, hdf5_err) - - ! Read global tallies - f_ptr = c_loc(global_tallies(1)) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - - ! Close global tallies dataset - call h5dclose_f(dset, hdf5_err) - - TALLIES_LOOP: do i = 1, n_tallies - ! Open tally group - call h5gopen_f(hdf5_state_point, "tallies/tally" // & - to_str(i), tally_group, hdf5_err) - - ! Read number of realizations - call hdf5_read_integer(tally_group, "n_realizations", & - tallies(i) % n_realizations) - - ! Open dataset for tally results - call h5dopen_f(tally_group, "results", dset, hdf5_err) - - ! Read sum and sum_sq for each tally bin - f_ptr = c_loc(tallies(i) % results(1,1)) - call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - - ! Close dataset for tally results - call h5dclose_f(dset, hdf5_err) - - ! Close tally group - call h5gclose_f(tally_group, hdf5_err) - end do TALLIES_LOOP - end if - - ! TODO: Use parallel HDF5 to read source bank in parallel - - ! Open dataset for source bank - call h5dopen_f(hdf5_state_point, "source_bank", dset, hdf5_err) - - ! Read source bank - f_ptr = c_loc(source_bank(1)) - call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - - ! Close dataset for source bank - call h5dclose_f(dset, hdf5_err) - - ! Close HDF5 state point file - call h5fclose_f(hdf5_state_point, hdf5_err) - - end subroutine hdf5_load_state_point + end subroutine hdf5_close_group !=============================================================================== -! HDF5_WRITE_INTEGER +! HDF5_WRITE_INTEGER writes integer scalar data !=============================================================================== subroutine hdf5_write_integer(group, name, buffer) - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer, intent(in) :: buffer + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: buffer ! data to write - integer :: rank = 1 - integer(HSIZE_T) :: dims(1) = (/1/) + ! Set rank and dimensions + hdf5_rank = 1 + dims1(1) = 1 - call h5ltmake_dataset_int_f(group, name, rank, dims, & + call h5ltmake_dataset_int_f(group, name, hdf5_rank, dims1, & (/ buffer /), hdf5_err) end subroutine hdf5_write_integer !=============================================================================== -! HDF5_WRITE_LONG +! HDF5_WRITE_INTEGER_1DARRAY writes integer 1-D array !=============================================================================== - subroutine hdf5_write_long(group, name, buffer) + subroutine hdf5_write_integer_1Darray(group, name, buffer, len) - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer(8), target, intent(in) :: buffer + integer, intent(in) :: len ! length of array to write + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: buffer(:) ! data to write - integer :: rank = 1 - integer(HSIZE_T) :: dims(1) = (/1/) - integer(HID_T) :: dspace - integer(HID_T) :: dset - type(c_ptr) :: f_ptr + ! Set rank and dimensions of data + hdf5_rank = 1 + dims1(1) = len + + ! Write data + call h5ltmake_dataset_int_f(group, name, hdf5_rank, dims1, & + buffer, hdf5_err) + + end subroutine hdf5_write_integer_1Darray + +!=============================================================================== +! HDF5_WRITE_INTEGER_2DARRAY write integer 2-D array +!=============================================================================== + + subroutine hdf5_write_integer_2Darray(group, name, buffer, length) + + integer, intent(in) :: length(2) ! length of array dimensions + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: buffer(length(1),length(2)) ! data to write + + ! Set rank and dimensions + hdf5_rank = 2 + dims2 = length + + ! Write data + call h5ltmake_dataset_int_f(group, name, hdf5_rank, dims2, & + buffer, hdf5_err) + + end subroutine hdf5_write_integer_2Darray + +!=============================================================================== +! HDF5_WRITE_INTEGER_3DARRAY writes integer 3-D array +!=============================================================================== + + subroutine hdf5_write_integer_3Darray(group, name, buffer, length) + + integer, intent(in) :: length(3) ! length of array dimensions + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: buffer(length(1),length(2), length(3)) ! data + + ! Set rank and dimensions + hdf5_rank = 3 + dims3 = length + + ! Write data + call h5ltmake_dataset_int_f(group, name, hdf5_rank, dims3, & + buffer, hdf5_err) + + end subroutine hdf5_write_integer_3Darray + +!=============================================================================== +! HDF5_WRITE_LONG writes long integer scalar data +!=============================================================================== + + subroutine hdf5_write_long(group, name, buffer, long_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer(8), target, intent(in) :: buffer ! data to write + integer(HID_T), intent(in) :: long_type ! HDF5 long type + + ! Set up rank and dimensions + hdf5_rank = 1 + dims1(1) = 1 ! Create dataspace and dataset - call h5screate_simple_f(rank, dims, dspace, hdf5_err) - call h5dcreate_f(group, name, hdf5_integer8_t, dspace, dset, hdf5_err) + call h5screate_simple_f(hdf5_rank, dims1, dspace, hdf5_err) + call h5dcreate_f(group, name, long_type, dspace, dset, hdf5_err) ! Write eight-byte integer f_ptr = c_loc(buffer) - call h5dwrite_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) + call h5dwrite_f(dset, long_type, f_ptr, hdf5_err) ! Close dataspace and dataset for long integer call h5dclose_f(dset, hdf5_err) @@ -1309,53 +299,204 @@ contains end subroutine hdf5_write_long !=============================================================================== -! HDF5_WRITE_DOUBLE +! HDF5_WRITE_DOUBLE writes double precision scalar data !=============================================================================== subroutine hdf5_write_double(group, name, buffer) - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - real(8), intent(in) :: buffer + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + real(8), intent(in) :: buffer ! data to write - integer :: rank = 1 - integer(HSIZE_T) :: dims(1) = (/1/) + ! Set up rank and dimensions + hdf5_rank = 1 + dims1(1) = 1 - call h5ltmake_dataset_double_f(group, name, rank, dims, & + ! Write data + call h5ltmake_dataset_double_f(group, name, hdf5_rank, dims1, & (/ buffer /), hdf5_err) end subroutine hdf5_write_double !=============================================================================== -! HDF5_READ_INTEGER +! HDF5_WRITE_DOUBLE_1DARRAY writes double precision 1-D array +!=============================================================================== + + subroutine hdf5_write_double_1Darray(group, name, buffer, len) + + integer, intent(in) :: len ! length of array + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + real(8), intent(in) :: buffer(:) ! data to write + + ! Set rank and dimensions of data + hdf5_rank = 1 + dims1(1) = len + + ! Write data + call h5ltmake_dataset_double_f(group, name, hdf5_rank, dims1, & + buffer, hdf5_err) + + end subroutine hdf5_write_double_1Darray + +!=============================================================================== +! HDF5_WRITE_DOUBLE_2DARRAY writes double precision 2-D array +!=============================================================================== + + subroutine hdf5_write_double_2Darray(group, name, buffer, length) + + integer, intent(in) :: length(2) ! length of array dimensions + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + real(8), intent(in) :: buffer(length(1),length(2)) ! data to write + + ! Set rank and dimensions of data + hdf5_rank = 2 + dims2 = length + + ! Write data + call h5ltmake_dataset_double_f(group, name, hdf5_rank, dims2, & + buffer, hdf5_err) + + end subroutine hdf5_write_double_2Darray + +!=============================================================================== +! HDF5_WRITE_DOUBLE_3DARRAY writes double precision 3-D aray +!=============================================================================== + + subroutine hdf5_write_double_3Darray(group, name, buffer, length) + + integer, intent(in) :: length(3) ! length of array dimensions + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + real(8), intent(in) :: buffer(length(1),length(2), length(3)) ! data + + ! Set rank and dimensions + hdf5_rank = 3 + dims3 = length + + ! Write data + call h5ltmake_dataset_double_f(group, name, hdf5_rank, dims3, & + buffer, hdf5_err) + + end subroutine hdf5_write_double_3Darray + +!=============================================================================== +! HDF5_WRITE_STRING writes string data +!=============================================================================== + + subroutine hdf5_write_string(group, name, buffer, length) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + character(*), intent(in) :: buffer ! data to write + integer, intent(in) :: length + + character(len=length), dimension(1) :: str_tmp + +! Fortran 2003 implementation not compatible with IBM compiler Feb 2013 +! type(c_ptr), dimension(1), target :: wdata +! character(len=length, kind=c_char), dimension(1), target :: c_str +! dims1(1) = 1 +! call h5screate_simple_f(1, dims1, dspace, hdf5_err) +! call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) +! c_str(1) = buffer +! wdata(1) = c_loc(c_str(1)) +! f_ptr = c_loc(wdata(1)) + + ! Number of strings to write + dims1(1) = 1 + + ! Insert null character at end of string when writing + call h5tset_strpad_f(H5T_STRING, H5T_STR_NULLPAD_F, hdf5_err) + + ! Create the dataspace and dataset + call h5screate_simple_f(1, dims1, dspace, hdf5_err) + call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) + + ! Set up dimesnions of string to write + dims2 = (/length, 1/) ! full array of strings to write + dims1(1) = length ! length of string + + ! Copy over string buffer to a rank 1 array + str_tmp(1) = buffer + + ! Write the variable dataset + call h5dwrite_vl_f(dset, H5T_STRING, str_tmp, dims2, dims1, hdf5_err, & + mem_space_id=dspace) + + ! Close all + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + + end subroutine hdf5_write_string + +!=============================================================================== +! HDF5_WRITE_ATTRIBUTE_STRING writes a string attribute to a variables +!=============================================================================== + + subroutine hdf5_write_attribute_string(group, var, attr_type, attr_str) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: var ! name of varaible to set attr + character(*), intent(in) :: attr_type ! the attr type id + character(*), intent(in) :: attr_str ! attribute sting + + call h5ltset_attribute_string_f(group, var, attr_type, attr_str, hdf5_err) + + end subroutine hdf5_write_attribute_string + +!=============================================================================== +! HDF5_READ_INTEGER reads integer scalar data !=============================================================================== subroutine hdf5_read_integer(group, name, buffer) - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer, intent(inout) :: buffer + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(inout) :: buffer ! read data to here - integer :: buffer_copy(1) - integer(HSIZE_T) :: dims(1) = (/1/) + integer :: buffer_copy(1) ! need an array for read - call h5ltread_dataset_int_f(group, name, buffer_copy, dims, hdf5_err) + ! Set up dimensions + dims1(1) = 1 + + ! Read data + call h5ltread_dataset_int_f(group, name, buffer_copy, dims1, hdf5_err) buffer = buffer_copy(1) end subroutine hdf5_read_integer !=============================================================================== -! HDF5_READ_LONG +! HDF5_READ_INTEGER_1DARRAY reads integer 1-D array !=============================================================================== - subroutine hdf5_read_long(group, name, buffer) + subroutine hdf5_read_integer_1Darray(group, name, buffer, length) - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - integer(8), target, intent(out) :: buffer + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(inout) :: buffer(:) ! read data to here + integer, intent(in) :: length ! length of array - integer(HID_T) :: dset - type(c_ptr) :: f_ptr + ! Set dimensions + dims1(1) = length + + ! Read data + call h5ltread_dataset_int_f(group, name, buffer, dims1, hdf5_err) + + end subroutine hdf5_read_integer_1Darray + + +!=============================================================================== +! HDF5_READ_LONG read long integer scalar data +!=============================================================================== + + subroutine hdf5_read_long(group, name, buffer, long_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer(8), target, intent(out) :: buffer ! read data to here + integer(HID_T), intent(in) :: long_type ! long integer type ! Open dataset call h5dopen_f(group, name, dset, hdf5_err) @@ -1364,7 +505,7 @@ contains f_ptr = c_loc(buffer) ! Read data from dataset - call h5dread_f(dset, hdf5_integer8_t, f_ptr, hdf5_err) + call h5dread_f(dset, long_type, f_ptr, hdf5_err) ! Close dataset call h5dclose_f(dset, hdf5_err) @@ -1372,23 +513,272 @@ contains end subroutine hdf5_read_long !=============================================================================== -! HDF5_READ_DOUBLE +! HDF5_READ_DOUBLE reads double precision scalar data !=============================================================================== subroutine hdf5_read_double(group, name, buffer) - integer(HID_T), intent(in) :: group - character(*), intent(in) :: name - real(8), intent(out) :: buffer + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + real(8), intent(inout) :: buffer ! read data to here - real(8) :: buffer_copy(1) - integer(HSIZE_T) :: dims(1) = (/1/) + real(8) :: buffer_copy(1) ! need an array for read - call h5ltread_dataset_double_f(group, name, buffer_copy, dims, hdf5_err) + ! Set up dimensions + dims1(1) = 1 + + ! Read data + call h5ltread_dataset_double_f(group, name, buffer_copy, dims1, hdf5_err) buffer = buffer_copy(1) end subroutine hdf5_read_double +!=============================================================================== +! HDF5_READ_DOUBLE_1DARRAY reads double precision 1-D array +!=============================================================================== + + subroutine hdf5_read_double_1Darray(group, name, buffer, length) + + integer(HID_T), intent(in) :: group ! name of group + integer, intent(in) :: length ! length of array + character(*), intent(in) :: name ! name of data + real(8), intent(inout) :: buffer(:) ! read data to here + + ! Set dimensions of data + dims1(1) = length + + ! Read data + call h5ltread_dataset_double_f(group, name, buffer, dims1, hdf5_err) + + end subroutine hdf5_read_double_1Darray + +!=============================================================================== +! HDF5_READ_STRING reads string data +!=============================================================================== + + subroutine hdf5_read_string(group, name, buffer) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + character(*), intent(inout) :: buffer ! read data to here + + call h5ltread_dataset_string_f(group, name, buffer, hdf5_err) + + end subroutine hdf5_read_string + +# ifdef MPI +!=============================================================================== +! HDF5_PARALLEL_READ_INTEGER reads interger scalar data in parallel +!=============================================================================== + + subroutine hdf5_parallel_read_integer(group, name, buffer, p_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, target, intent(inout) :: buffer ! read data to here + integer, intent(in) :: p_type ! independent or collective I/O + + ! Set up dimension + dims1(1) = 1 + + ! Create property list for independent or collective read + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + + ! Set independent or collective option + call h5pset_dxpl_mpio_f(plist, p_type, hdf5_err) + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Read data + f_ptr = c_loc(buffer) + call h5dread_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) + + ! Close dataset and property list + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_read_integer + +!=============================================================================== +! HDF5_PARALLEL_READ_INTEGER_1DARRAY reads integer 1-D array in parallel +!=============================================================================== + + subroutine hdf5_parallel_read_integer_1Darray(group, name, buffer, length, & + p_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: length ! length of array + integer, target, intent(inout) :: buffer(length) ! read data to here + integer, intent(in) :: p_type ! independent or collective I/O + + ! Create property list for independent or collective read + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + + ! Set independent or collective option + call h5pset_dxpl_mpio_f(plist, p_type, hdf5_err) + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Read data + f_ptr = c_loc(buffer(1)) + call h5dread_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, xfer_prp=plist) + + ! Close dataset and property list + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_read_integer_1Darray + + +!=============================================================================== +! HDF5_PARALLEL_READ_LONG read long integer scalar data in parallel +!=============================================================================== + + subroutine hdf5_parallel_read_long(group, name, buffer, long_type, p_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer(8), target, intent(out) :: buffer ! read data to here + integer(HID_T), intent(in) :: long_type ! long integer type + integer, intent(in) :: p_type ! independent or collective I/O + + ! Create property list for independent or collective read + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + + ! Set independent or collective option + call h5pset_dxpl_mpio_f(plist, p_type, hdf5_err) + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Read data + f_ptr = c_loc(buffer) + call h5dread_f(dset, long_type, f_ptr, hdf5_err, xfer_prp=plist) + + ! Close dataset and property list + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_read_long + +!=============================================================================== +! HDF5_PARALLEL_READ_DOUBLE reads double precision scalar data in parallel +!=============================================================================== + + subroutine hdf5_parallel_read_double(group, name, buffer, p_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + real(8), target, intent(inout) :: buffer ! read data to here + integer, intent(in) :: p_type ! independent or collective I/O + + ! Create property list for independent or collective read + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + + ! Set independent or collective option + call h5pset_dxpl_mpio_f(plist, p_type, hdf5_err) + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Read data + f_ptr = c_loc(buffer) + call h5dread_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) + + ! Close dataset and property list + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_read_double + +!=============================================================================== +! HDF5_PARLLEL_READ_DOUBLE_1DARRAY reads double precision 1-D array in parallel +!=============================================================================== + + subroutine hdf5_parallel_read_double_1Darray(group, name, buffer, length, & + p_type) + + integer(HID_T), intent(in) :: group ! name of group + integer, intent(in) :: length ! length of array + integer, intent(in) :: p_type ! indepedent or collective I/O + character(*), intent(in) :: name ! name of data + real(8), target, intent(inout) :: buffer(length) ! read data to here + + ! Create property list for independent or collective read + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + + ! Set independent or collective option + call h5pset_dxpl_mpio_f(plist, p_type, hdf5_err) + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Read data + f_ptr = c_loc(buffer(1)) + call h5dread_f(dset, H5T_NATIVE_DOUBLE, f_ptr, hdf5_err, xfer_prp=plist) + + ! Close dataset and property list + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_read_double_1Darray + +!=============================================================================== +! HDF5_PARALLEL_READ_STRING reads string data +!=============================================================================== + + subroutine hdf5_parallel_read_string(group, name, buffer, length, p_type) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + character(*), target, intent(inout) :: buffer ! read data to here + integer, intent(in) :: p_type ! independent or collective IO + integer, intent(in) :: length ! length of string + + character(len=length), dimension(1) :: str_tmp + ! Fortran 2003 implementation not compatible with IBM Feb 2013 compiler +! type(c_ptr), dimension(1), target :: buf_ptr +! character(len=length, kind=c_char), pointer :: chr_ptr +! f_ptr = c_loc(buf_ptr(1)) +! call h5dread_f(dset, H5T_STRING, f_ptr, hdf5_err, xfer_prp=plist) +! call c_f_pointer(buf_ptr(1), chr_ptr) +! buffer = chr_ptr +! nullify(chr_ptr) + + ! Create property list for independent or collective read + call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err) + + ! Set independent or collective option + call h5pset_dxpl_mpio_f(plist, p_type, hdf5_err) + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Get dataspace to read + call h5dget_space_f(dset, dspace, hdf5_err) + + ! Set dimensions + dims2 = (/length, 1/) + dims1(1) = length + + ! Read in the data + call h5dread_vl_f(dset, H5T_STRING, str_tmp, dims2, dims1, hdf5_err, & + mem_space_id=dspace, xfer_prp = plist) + + ! Copy over buffer + buffer = str_tmp(1) + + ! Close dataset and property list + call h5dclose_f(dset, hdf5_err) + call h5pclose_f(plist, hdf5_err) + + end subroutine hdf5_parallel_read_string + +# endif + #endif end module hdf5_interface diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 new file mode 100644 index 0000000000..3963b0035b --- /dev/null +++ b/src/hdf5_summary.F90 @@ -0,0 +1,715 @@ +module hdf5_summary + +#ifdef HDF5 + + use ace_header, only: Reaction, UrrData, Nuclide + use constants + use endf, only: reaction_name + use geometry_header, only: Cell, Surface, Universe, Lattice + use global + use material_header, only: Material + use mesh_header, only: StructuredMesh + use output_interface + use output, only: time_stamp + use string, only: to_str + use tally_header, only: TallyObject + +contains + +!=============================================================================== +! HDF5_WRITE_SUMMARY +!=============================================================================== + + subroutine hdf5_write_summary() + + character(MAX_FILE_LEN) :: filename = "summary.h5" + + ! Create a new file using default properties. + call file_create(filename, "serial") + + ! Write header information + call hdf5_write_header() + + ! Write eigenvalue information + if (run_mode == MODE_EIGENVALUE) then + + ! Write number of particles + call write_data(n_particles, "n_particles") + + ! Use H5LT interface to write n_batches, n_inactive, and n_active + call write_data(n_batches, "n_batches") + call write_data(n_inactive, "n_inactive") + call write_data(n_active, "n_active") + call write_data(gen_per_batch, "gen_per_batch") + + ! Add description of each variable + call write_attribute_string("n_particles", & + "description", "Number of particles per generation") + call write_attribute_string("n_batches", & + "description", "Total number of batches") + call write_attribute_string("n_inactive", & + "description", "Number of inactive batches") + call write_attribute_string("n_active", & + "description", "Number of active batches") + call write_attribute_string("gen_per_batch", & + "description", "Number of generations per batch") + end if + + call hdf5_write_geometry() + call hdf5_write_materials() + call hdf5_write_nuclides() + if (n_tallies > 0) then + call hdf5_write_tallies() + end if + + ! Terminate access to the file. + call file_close("serial") + + end subroutine hdf5_write_summary + +!=============================================================================== +! HDF5_WRITE_HEADER +!=============================================================================== + + subroutine hdf5_write_header() + + ! Write version information + call write_data(VERSION_MAJOR, "version_major") + call write_data(VERSION_MINOR, "version_minor") + call write_data(VERSION_RELEASE, "version_release") + + ! Write current date and time + call write_data(time_stamp(), "date_and_time") + + ! Write MPI information + call write_data(n_procs, "n_procs") + call write_attribute_string("n_procs", "description", & + "Number of MPI processes") + + end subroutine hdf5_write_header + +!=============================================================================== +! HDF5_WRITE_GEOMETRY +!=============================================================================== + + subroutine hdf5_write_geometry() + + integer :: i, j, k, m + integer :: n_x, n_y, n_z + integer, allocatable :: lattice_universes(:,:,:) + type(Cell), pointer :: c => null() + type(Surface), pointer :: s => null() + type(Universe), pointer :: u => null() + type(Lattice), pointer :: lat => null() + + ! Use H5LT interface to write number of geometry objects + call write_data(n_cells, "n_cells", group="geometry") + call write_data(n_surfaces, "n_surfaces", group="geometry") + call write_data(n_universes, "n_universes", group="geometry") + call write_data(n_lattices, "n_lattices", group="geometry") + + ! ========================================================================== + ! WRITE INFORMATION ON CELLS + + ! Create a cell group (nothing directly written in this group) then close + call hdf5_open_group("geometry/cells") + call hdf5_close_group() + + ! Write information on each cell + CELL_LOOP: do i = 1, n_cells + c => cells(i) + + ! Write universe for this cell + call write_data(universes(c % universe) % id, "universe", & + group="geometry/cells/cell " // trim(to_str(c % id))) + + ! Write information on what fills this cell + select case (c % type) + case (CELL_NORMAL) + call write_data("normal", "fill_type", & + group="geometry/cells/cell " // trim(to_str(c % id))) + if (c % material == MATERIAL_VOID) then + call write_data(-1, "material", & + group="geometry/cells/cell " // trim(to_str(c % id))) + else + call write_data(materials(c % material) % id, "material", & + group="geometry/cells/cell " // trim(to_str(c % id))) + end if + case (CELL_FILL) + call write_data("universe", "fill_type", & + group="geometry/cells/cell " // trim(to_str(c % id))) + call write_data(universes(c % fill) % id, "material", & + group="geometry/cells/cell " // trim(to_str(c % id))) + case (CELL_LATTICE) + call write_data("lattice", "fill_type", & + group="geometry/cells/cell " // trim(to_str(c % id))) + call write_data(lattices(c % fill) % id, "lattice", & + group="geometry/cells/cell " // trim(to_str(c % id))) + end select + + ! Write list of bounding surfaces + if (c % n_surfaces > 0) then + call write_data(c % surfaces, "surfaces", length= c % n_surfaces, & + group="geometry/cells/cell " // trim(to_str(c % id))) + end if + + end do CELL_LOOP + + ! ========================================================================== + ! WRITE INFORMATION ON SURFACES + + ! Create surfaces group (nothing directly written here) then close + call hdf5_open_group("geometry/surfaces") + call hdf5_close_group() + + ! Write information on each surface + SURFACE_LOOP: do i = 1, n_surfaces + s => surfaces(i) + + ! Write surface type + select case (s % type) + case (SURF_PX) + call write_data("X Plane", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_PY) + call write_data("Y Plane", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_PZ) + call write_data("Z Plane", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_PLANE) + call write_data("Plane", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_CYL_X) + call write_data("X Cylinder", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_CYL_Y) + call write_data("Y Cylinder", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_CYL_Z) + call write_data("Z Cylinder", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_SPHERE) + call write_data("Sphere", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_CONE_X) + call write_data("X Cone", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_CONE_Y) + call write_data("Y Cone", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (SURF_CONE_Z) + call write_data("Z Cone", "type", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + end select + + ! Write coefficients for surface + call write_data(s % coeffs, "coefficients", length=size(s % coeffs), & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + + ! Write positive neighbors + if (allocated(s % neighbor_pos)) then + call write_data(s % neighbor_pos, "neighbors_positive", & + length=size(s % neighbor_pos), & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + end if + + ! Write negative neighbors + if (allocated(s % neighbor_neg)) then + call write_data(s % neighbor_neg, "neighbors_negative", & + length=size(s % neighbor_neg), & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + end if + + ! Write boundary condition + select case (s % bc) + case (BC_TRANSMIT) + call write_data("transmission", "boundary_condition", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (BC_VACUUM) + call write_data("vacuum", "boundary_condition", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (BC_REFLECT) + call write_data("reflective", "boundary_condition", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + case (BC_PERIODIC) + call write_data("periodic", "boundary_condition", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + end select + + end do SURFACE_LOOP + + ! ========================================================================== + ! WRITE INFORMATION ON UNIVERSES + + ! Create universes group (nothing directly written here) then close + call hdf5_open_group("geometry/universes") + call hdf5_close_group() + + ! Write information on each universe + UNIVERSE_LOOP: do i = 1, n_universes + u => universes(i) + + ! Write list of cells in this universe + if (u % n_cells > 0) then + call write_data(u % cells, "cells", length=u % n_cells, & + group="geometry/universes/universe " // trim(to_str(u % id))) + end if + + end do UNIVERSE_LOOP + + ! ========================================================================== + ! WRITE INFORMATION ON LATTICES + + ! Create lattices group (nothing directly written here) then close + call hdf5_open_group("geometry/lattices") + call hdf5_close_group() + + ! Write information on each lattice + LATTICE_LOOP: do i = 1, n_lattices + lat => lattices(i) + + ! Write lattice type + select case(lat % type) + case (LATTICE_RECT) + call write_data("rectangular", "type", & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + case (LATTICE_HEX) + call write_data("hexagonal", "type", & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + end select + + ! Write lattice dimensions, lower left corner, and width of element + call write_data(lat % dimension, "dimension", & + length=lat % n_dimension, & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_data(lat % lower_left, "lower_left", & + length=lat % n_dimension, & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call write_data(lat % width, "width", & + length=lat % n_dimension, & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + + ! Determine dimensions of lattice + n_x = lat % dimension(1) + n_y = lat % dimension(2) + if (lat % n_dimension == 3) then + n_z = lat % dimension(3) + else + n_z = 1 + end if + + ! Write lattice universes + allocate(lattice_universes(n_x, n_y, n_z)) + do j = 1, n_x + do k = 1, n_y + do m = 1, n_z + lattice_universes(j,k,m) = universes(lat % universes(j,k,m)) % id + end do + end do + end do + call write_data(lattice_universes, "universes", & + length=(/n_x, n_y, n_z/), & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + deallocate(lattice_universes) + + end do LATTICE_LOOP + + end subroutine hdf5_write_geometry + +!=============================================================================== +! HDF5_WRITE_MATERIALS +!=============================================================================== + + subroutine hdf5_write_materials() + + integer :: i + integer :: j + integer, allocatable :: zaids(:) + type(Material), pointer :: m => null() + + ! Use H5LT interface to write number of materials + call write_data(n_materials, "n_materials", group="materials") + + ! Write information on each material + do i = 1, n_materials + m => materials(i) + + ! Write atom density with units + call write_data(m % density, "atom_density", & + group="materials/material " // trim(to_str(m % id))) + call write_attribute_string("atom_density", "units", "atom/b-cm", & + group="materials/material " // trim(to_str(m % id))) + + ! Copy ZAID for each nuclide to temporary array + allocate(zaids(m % n_nuclides)) + do j = 1, m % n_nuclides + zaids(j) = nuclides(m % nuclide(j)) % zaid + end do + + ! Write temporary array to 'nuclides' + call write_data(zaids, "nuclides", length=m % n_nuclides, & + group="materials/material " // trim(to_str(m % id))) + + ! Deallocate temporary array + deallocate(zaids) + + ! Write atom densities + call write_data(m % atom_density, "nuclide_densities", & + length=m % n_nuclides, & + group="materials/material " // trim(to_str(m % id))) + + ! Write S(a,b) information if present + if (m % n_sab > 0) then + call write_data(m % i_sab_nuclides, "i_sab_nuclides", & + length=m % n_sab, & + group="materials/material " // trim(to_str(m % id))) + call write_data(m % i_sab_tables, "i_sab_tables", & + length=m % n_sab, & + group="materials/material " // trim(to_str(m % id))) + end if + + end do + + end subroutine hdf5_write_materials + +!=============================================================================== +! HDF5_WRITE_TALLIES +!=============================================================================== + + subroutine hdf5_write_tallies() + + integer :: i, j + integer, allocatable :: temp_array(:) ! nuclide bin array + type(StructuredMesh), pointer :: m => null() + type(TallyObject), pointer :: t => null() + + ! Write total number of meshes + call write_data(n_meshes, "n_meshes", group="tallies") + + ! Write information for meshes + MESH_LOOP: do i = 1, n_meshes + m => meshes(i) + + ! Write type and number of dimensions + call write_data(m % type, "type", & + group="tallies/mesh " // trim(to_str(m % id))) + + call write_data(m % n_dimension, "n_dimension", & + group="tallies/mesh " // trim(to_str(m % id))) + + ! Write mesh information + call write_data(m % dimension, "dimension", & + length=m % n_dimension, & + group="tallies/mesh " // trim(to_str(m % id))) + call write_data(m % lower_left, "lower_left", & + length=m % n_dimension, & + group="tallies/mesh " // trim(to_str(m % id))) + call write_data(m % upper_right, "upper_right", & + length=m % n_dimension, & + group="tallies/mesh " // trim(to_str(m % id))) + call write_data(m % width, "width", & + length=m % n_dimension, & + group="tallies/mesh " // trim(to_str(m % id))) + + end do MESH_LOOP + + ! Write number of tallies + call write_data(n_tallies, "n_tallies", group="tallies") + + TALLY_METADATA: do i = 1, n_tallies + ! Get pointer to tally + t => tallies(i) + + ! Write size of each tally + call write_data(t % total_score_bins, "total_score_bins", & + group="tallies/tally " // trim(to_str(t % id))) + call write_data(t % total_filter_bins, "total_filter_bins", & + group="tallies/tally " // trim(to_str(t % id))) + + ! Write number of filters + call write_data(t % n_filters, "n_filters", & + group="tallies/tally " // trim(to_str(t % id))) + + FILTER_LOOP: do j = 1, t % n_filters + ! Write type of filter + call write_data(t % filters(j) % type, "type", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + + ! Write number of bins for this filter + call write_data(t % filters(j) % n_bins, "n_bins", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + + ! Write filter bins + if (t % filters(j) % type == FILTER_ENERGYIN .or. & + t % filters(j) % type == FILTER_ENERGYOUT) then + call write_data(t % filters(j) % real_bins, "bins", & + length=size(t % filters(j) % real_bins), & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + else + call write_data(t % filters(j) % int_bins, "bins", & + length=size(t % filters(j) % int_bins), & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + end if + + ! Write name of type + select case (t % filters(j) % type) + case(FILTER_UNIVERSE) + call write_data("universe", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_MATERIAL) + call write_data("material", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_CELL) + call write_data("cell", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_CELLBORN) + call write_data("cellborn", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_SURFACE) + call write_data("surface", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_MESH) + call write_data("mesh", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_ENERGYIN) + call write_data("energy", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + case(FILTER_ENERGYOUT) + call write_data("energyout", "type_name", & + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) + end select + + end do FILTER_LOOP + + ! Write number of nuclide bins + call write_data(t % n_nuclide_bins, "n_nuclide_bins", & + group="tallies/tally " // trim(to_str(t % id))) + + ! Create temporary array for nuclide bins + allocate(temp_array(t % n_nuclide_bins)) + NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins + if (t % nuclide_bins(j) > 0) then + temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid + else + temp_array(j) = t % nuclide_bins(j) + end if + end do NUCLIDE_LOOP + + ! Write and deallocate nuclide bins + call write_data(temp_array, "nuclide_bins", length=t % n_nuclide_bins, & + group="tallies/tally " // trim(to_str(t % id))) + deallocate(temp_array) + + ! Write number of score bins + call write_data(t % n_score_bins, "n_score_bins", & + group="tallies/tally " // trim(to_str(t % id))) + call write_data(t % score_bins, "score_bins", length=t % n_score_bins, & + group="tallies/tally " // trim(to_str(t % id))) + + end do TALLY_METADATA + + end subroutine hdf5_write_tallies + +!=============================================================================== +! HDF5_WRITE_NUCLIDES +!=============================================================================== + + subroutine hdf5_write_nuclides() + + integer :: i, j + integer :: size_total + integer :: size_xs + integer :: size_angle + integer :: size_energy + type(Nuclide), pointer :: nuc => null() + type(Reaction), pointer :: rxn => null() + type(UrrData), pointer :: urr => null() + + ! Use H5LT interface to write number of nuclides + call write_data(n_nuclides_total, "n_nuclides", group="nuclides") + + ! Write information on each nuclide + NUCLIDE_LOOP: do i = 1, n_nuclides_total + nuc => nuclides(i) + + ! Determine size of cross-sections + size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8 + size_total = size_xs + + ! Write some basic attributes + call write_data(nuc % zaid, "zaid", & + group="nuclides/" // trim(nuc % name)) + call write_data(nuc % awr, "awr", & + group="nuclides/" // trim(nuc % name)) + call write_data(nuc % kT, "kT", & + group="nuclides/" // trim(nuc % name)) + call write_data(nuc % n_grid, "n_grid", & + group="nuclides/" // trim(nuc % name)) + call write_data(nuc % n_reaction, "n_reactions", & + group="nuclides/" // trim(nuc % name)) + call write_data(nuc % n_fission, "n_fission", & + group="nuclides/" // trim(nuc % name)) + call write_data(size_xs, "size_xs", & + group="nuclides/" // trim(nuc % name)) + + ! ======================================================================= + ! WRITE INFORMATION ON EACH REACTION + + ! Create overall group for reactions and close it + call hdf5_open_group("nuclides/" // trim(nuc % name) // "/reactions") + call hdf5_close_group() + + RXN_LOOP: do j = 1, nuc % n_reaction + ! Information on each reaction + rxn => nuc % reactions(j) + + ! Determine size of angle distribution + if (rxn % has_angle_dist) then + size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 + else + size_angle = 0 + end if + + ! Determine size of energy distribution + if (rxn % has_energy_dist) then + size_energy = size(rxn % edist % data) * 8 + else + size_energy = 0 + end if + + ! Write information on reaction + call write_data(rxn % Q_value, "Q_value", & + group="nuclides/" // trim(nuc % name) // "/reactions/" // & + trim(reaction_name(rxn % MT))) + call write_data(rxn % multiplicity, "multiplicity", & + group="nuclides/" // trim(nuc % name) // "/reactions/" // & + trim(reaction_name(rxn % MT))) + call write_data(rxn % threshold, "threshold", & + group="nuclides/" // trim(nuc % name) // "/reactions/" // & + trim(reaction_name(rxn % MT))) + call write_data(size_angle, "size_angle", & + group="nuclides/" // trim(nuc % name) // "/reactions/" // & + trim(reaction_name(rxn % MT))) + call write_data(size_energy, "size_energy", & + group="nuclides/" // trim(nuc % name) // "/reactions/" // & + trim(reaction_name(rxn % MT))) + + ! Accumulate data size + size_total = size_total + size_angle + size_energy + end do RXN_LOOP + + ! ======================================================================= + ! WRITE INFORMATION ON URR PROBABILITY TABLES + + if (nuc % urr_present) then + urr => nuc % urr_data + call write_data(urr % n_energy, "urr_n_energy", & + group="nuclides/" // trim(nuc % name)) + call write_data(urr % n_prob, "urr_n_prob", & + group="nuclides/" // trim(nuc % name)) + call write_data(urr % interp, "urr_interp", & + group="nuclides/" // trim(nuc % name)) + call write_data(urr % inelastic_flag, "urr_inelastic", & + group="nuclides/" // trim(nuc % name)) + call write_data(urr % absorption_flag, "urr_absorption", & + group="nuclides/" // trim(nuc % name)) + call write_data(urr % energy(1), "urr_min_E", & + group="nuclides/" // trim(nuc % name)) + call write_data(urr % energy(urr % n_energy), "urr_max_E", & + group="nuclides/" // trim(nuc % name)) + end if + + ! Write total memory used + call write_data(size_total, "size_total", & + group="nuclides/" // trim(nuc % name)) + + end do NUCLIDE_LOOP + + end subroutine hdf5_write_nuclides + +!=============================================================================== +! HDF5_WRITE_TIMING +!=============================================================================== + + subroutine hdf5_write_timing() + + integer(8) :: total_particles + real(8) :: speed + + ! Write timing data + call write_data(time_initialize % elapsed, "time_initialize", & + group="timing") + call write_data(time_read_xs % elapsed, "time_read_xs", & + group="timing") + call write_data(time_unionize % elapsed, "time_unionize", & + group="timing") + call write_data(time_transport % elapsed, "time_transport", & + group="timing") + call write_data(time_bank % elapsed, "time_bank", & + group="timing") + call write_data(time_bank_sample % elapsed, "time_bank_sample", & + group="timing") + call write_data(time_bank_sendrecv % elapsed, "time_bank_sendrecv", & + group="timing") + call write_data(time_tallies % elapsed, "time_tallies", & + group="timing") + call write_data(time_inactive % elapsed, "time_inactive", & + group="timing") + call write_data(time_active % elapsed, "time_active", & + group="timing") + call write_data(time_finalize % elapsed, "time_finalize", & + group="timing") + call write_data(time_total % elapsed, "time_total", & + group="timing") + + ! Add descriptions to timing data + call write_attribute_string("time_initialize", "description", & + "Total time elapsed for initialization (s)", group="timing") + call write_attribute_string("time_read_xs", "description", & + "Time reading cross-section libraries (s)", group="timing") + call write_attribute_string("time_unionize", "description", & + "Time unionizing energy grid (s)", group="timing") + call write_attribute_string("time_transport", "description", & + "Time in transport only (s)", group="timing") + call write_attribute_string("time_bank", "description", & + "Total time synchronizing fission bank (s)", group="timing") + call write_attribute_string("time_bank_sample", "description", & + "Time between generations sampling source sites (s)", group="timing") + call write_attribute_string("time_bank_sendrecv", "description", & + "Time between generations SEND/RECVing source sites (s)", & + group="timing") + call write_attribute_string("time_tallies", "description", & + "Time between batches accumulating tallies (s)", group="timing") + call write_attribute_string("time_inactive", "description", & + "Total time in inactive batches (s)", group="timing") + call write_attribute_string("time_active", "description", & + "Total time in active batches (s)", group="timing") + call write_attribute_string("time_finalize", "description", & + "Total time for finalization (s)", group="timing") + call write_attribute_string("time_total", "description", & + "Total time elapsed (s)", group="timing") + + ! Write calculation rate + total_particles = n_particles * n_batches * gen_per_batch + speed = real(total_particles) / (time_inactive % elapsed + & + time_active % elapsed) + call write_data(speed, "neutrons_per_second", group="timing") + + end subroutine hdf5_write_timing + +#endif + +end module hdf5_summary diff --git a/src/initialize.F90 b/src/initialize.F90 index dd28f9bf9c..a04c8a2200 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -13,11 +13,12 @@ module initialize cells_in_univ_dict, read_plots_xml use output, only: title, header, write_summary, print_version, & print_usage, write_xs_summary, print_plot + use output_interface, only: file_open, file_close, read_data use random_lcg, only: initialize_prng use source, only: initialize_source use state_point, only: load_state_point use string, only: to_str, str_to_int, starts_with, ends_with - use tally_header, only: TallyObject + use tally_header, only: TallyObject, TallyResult use tally_initialize, only: configure_tallies #ifdef MPI @@ -25,8 +26,8 @@ module initialize #endif #ifdef HDF5 - use hdf5_interface, only: hdf5_initialize, hdf5_write_summary, & - hdf5_load_state_point + use hdf5_interface + use hdf5_summary, only: hdf5_write_summary #endif implicit none @@ -121,11 +122,7 @@ contains ! If this is a restart run, load the state point data and binary source ! file -#ifdef HDF5 - if (restart_run) call hdf5_load_state_point() -#else if (restart_run) call load_state_point() -#endif end if if (master) then @@ -239,6 +236,52 @@ contains end subroutine initialize_mpi #endif +#ifdef HDF5 + +!=============================================================================== +! HDF5_INITIALIZE +!=============================================================================== + + subroutine hdf5_initialize() + + type(TallyResult), target :: tmp(2) ! temporary TallyResult + type(Bank), target :: tmpb(2) ! temporary Bank + 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 the compound datatype for TallyResult + call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), & + c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err) + call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), & + c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err) + call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), & + c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, 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) + + ! Determine type for integer(8) + hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND) + + end subroutine hdf5_initialize + +#endif + !=============================================================================== ! READ_COMMAND_LINE reads all parameters from the command line !=============================================================================== @@ -248,6 +291,7 @@ contains integer :: i ! loop index integer :: argc ! number of command line arguments integer :: last_flag ! index of last flag + integer :: filetype character(MAX_FILE_LEN) :: pwd ! present working directory character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments @@ -284,10 +328,23 @@ contains call fatal_error() end if case ('-r', '-restart', '--restart') - ! Read path for state point + ! Read path for state point/particle restart i = i + 1 - path_state_point = argv(i) - restart_run = .true. + + ! Check what type of file this is + call file_open(argv(i), 'parallel', 'r') + call read_data(filetype, 'filetype') + call file_close('parallel') + + ! Set path and flag for type of run + select case (filetype) + case (FILETYPE_STATEPOINT) + path_state_point = argv(i) + restart_run = .true. + case (FILETYPE_PARTICLE_RESTART) + path_particle_restart = argv(i) + particle_restart_run = .true. + end select case ('-t', '-tallies', '--tallies') run_mode = MODE_TALLIES @@ -306,11 +363,6 @@ contains case ('-eps_tol', '-ksp_gmres_restart') ! Handle options that would be based to PETSC i = i + 1 - case ('-s','-particle','--particle') - ! Read in path for particle restart - i = i + 1 - path_particle_restart = argv(i) - particle_restart_run = .true. case default message = "Unknown command line option: " // argv(i) call fatal_error() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 75d0518d1f..a9b62d4905 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -148,7 +148,7 @@ contains gen_per_batch = eigenvalue_ % generations_per_batch ! Allocate array for batch keff and entropy - allocate(k_batch(n_batches)) + allocate(k_generation(n_batches*gen_per_batch)) allocate(entropy(n_batches*gen_per_batch)) entropy = ZERO end if diff --git a/src/mpiio_interface.F90 b/src/mpiio_interface.F90 new file mode 100644 index 0000000000..d3e9774df8 --- /dev/null +++ b/src/mpiio_interface.F90 @@ -0,0 +1,239 @@ +module mpiio_interface + +#ifdef MPI + use mpi + + implicit none + + integer :: mpi_fh ! MPI file handle + integer :: mpiio_err ! MPI error code + +contains + +!=============================================================================== +! MPI_CREATE_FILE creates a file using MPI file I/O +!=============================================================================== + + subroutine mpi_create_file(filename, fh) + + character(*), intent(in) :: filename ! name of file to create + integer, intent(inout) :: fh ! file handle + + ! Create the file + call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, MPI_MODE_CREATE + & + MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpiio_err) + + end subroutine mpi_create_file + +!=============================================================================== +! MPI_OPEN_FILE opens a file using MPI file I/O +!=============================================================================== + + subroutine mpi_open_file(filename, fh, mode) + + character(*), intent(in) :: filename ! name of file to open + character(*), intent(in) :: mode ! open 'r' read, 'w' write + integer, intent(inout) :: fh ! file handle + + integer :: open_mode + + ! Determine access mode + open_mode = MPI_MODE_RDONLY + if (mode == 'w') then + open_mode = MPI_MODE_WRONLY + end if + + ! Create the file + call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, & + open_mode, MPI_INFO_NULL, fh, mpiio_err) + + end subroutine mpi_open_file + +!=============================================================================== +! MPI_CLOSE_FILE closes a file using MPI file I/O +!=============================================================================== + + subroutine mpi_close_file(fh) + + integer, intent(inout) :: fh ! file handle + + call MPI_FILE_CLOSE(fh, mpiio_err) + + end subroutine mpi_close_file + +!=============================================================================== +! MPI_WRITE_INTEGER writes integer scalar data using MPI File I/O +!=============================================================================== + + subroutine mpi_write_integer(fh, buffer) + + integer, intent(in) :: fh ! file handle + integer, intent(in) :: buffer ! data to write + + call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_write_integer + +!=============================================================================== +! MPI_WRITE_INTEGER_1DARRAY writes integer 1-D array data using MPI File I/O +!=============================================================================== + + subroutine mpi_write_integer_1Darray(fh, buffer, length) + + integer, intent(in) :: fh ! file handle + integer, intent(in) :: length ! length of array + integer, intent(in) :: buffer(:) ! data to write + + call MPI_FILE_WRITE(fh, buffer, length, MPI_INTEGER, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_write_integer_1Darray + +!=============================================================================== +! MPI_WRITE_LONG writes long integer scalar data using MPI file I/O +!=============================================================================== + + subroutine mpi_write_long(fh, buffer) + + integer, intent(in) :: fh ! file handle + integer(8), intent(in) :: buffer ! data to write + + call MPI_FILE_WRITE(fh, buffer, 1, MPI_INTEGER8, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_write_long + +!=============================================================================== +! MPI_WRITE_DOUBLE writes double precision scalar data using MPI file I/O +!=============================================================================== + + subroutine mpi_write_double(fh, buffer) + + integer, intent(in) :: fh ! file handle + real(8), intent(in) :: buffer ! data to write + + call MPI_FILE_WRITE(fh, buffer, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_write_double + +!=============================================================================== +! MPI_WRITE_DOUBLE_1DARRAY writes double precision 1-D array using MPI file I/O +!=============================================================================== + + subroutine mpi_write_double_1Darray(fh, buffer, length) + + integer, intent(in) :: fh ! file handle + integer, intent(in) :: length ! length of array + real(8), intent(in) :: buffer(:) ! data to write + + call MPI_FILE_WRITE(fh, buffer, length, MPI_REAL8, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_write_double_1Darray + +!=============================================================================== +! MPI_WRITE_STRING writes string data using MPI file I/O +!=============================================================================== + + subroutine mpi_write_string(fh, buffer, length) + + character(*), intent(in) :: buffer ! data to write + integer, intent(in) :: fh ! file handle + integer, intent(in) :: length ! length of data + + call MPI_FILE_WRITE(fh, buffer, length, MPI_CHARACTER, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_write_string + +!=============================================================================== +! MPI_READ_INTEGER reads integer scalar data using MPI file I/O +!=============================================================================== + + subroutine mpi_read_integer(fh, buffer) + + integer, intent(in) :: fh ! file handle + integer, intent(inout) :: buffer ! read data to here + + call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_read_integer + +!=============================================================================== +! MPI_READ_INTEGER_1DARRAY reads integer 1-D array using MPI file I/O +!=============================================================================== + + subroutine mpi_read_integer_1Darray(fh, buffer, length) + + integer, intent(in) :: fh ! file handle + integer, intent(in) :: length ! length of array + integer, intent(inout) :: buffer(:) ! read data to here + + call MPI_FILE_READ(fh, buffer, length, MPI_INTEGER, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_read_integer_1Darray + +!=============================================================================== +! MPI_READ_LONG reads long integer scalar data using MPI file I/O +!=============================================================================== + + subroutine mpi_read_long(fh, buffer) + + integer, intent(in) :: fh ! file handle + integer(8), intent(inout) :: buffer ! read data to here + + call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER8, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_read_long + +!=============================================================================== +! MPI_READ_DOUBLE reads double precision scalar data using MPI file I/O +!=============================================================================== + + subroutine mpi_read_double(fh, buffer) + + integer, intent(in) :: fh ! file handle + real(8), intent(inout) :: buffer ! read data to here + + call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_read_double + +!=============================================================================== +! MPI_READ_DOUBLE_1DARRAY reads double precision 1-D array using MPI file I/O +!=============================================================================== + + subroutine mpi_read_double_1Darray(fh, buffer, length) + + integer, intent(in) :: fh ! file handle + integer, intent(in) :: length ! length of array + real(8), intent(inout) :: buffer(:) ! read data to here + + call MPI_FILE_READ(fh, buffer, length, MPI_REAL8, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_read_double_1Darray + +!=============================================================================== +! MPI_READ_STRING reads string data using MPI file I/O +!=============================================================================== + + subroutine mpi_read_string(fh, buffer, length) + + character(*), intent(inout) :: buffer ! read data to here + integer, intent(in) :: fh ! file handle + integer, intent(in) :: length ! length of string + + call MPI_FILE_READ(fh, buffer, length, MPI_CHARACTER, & + MPI_STATUS_IGNORE, mpiio_err) + + end subroutine mpi_read_string + +#endif +end module mpiio_interface diff --git a/src/output.F90 b/src/output.F90 index 4ca455624b..b8302226d7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -166,8 +166,8 @@ contains write(OUTPUT_UNIT,*) write(OUTPUT_UNIT,*) 'Options:' write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode' - write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run' - write(OUTPUT_UNIT,*) ' -s, --particle Run a single particle history' + write(OUTPUT_UNIT,*) ' -r, --restart Restart a previous run from a state point' + write(OUTPUT_UNIT,*) ' or a particle restart file' write(OUTPUT_UNIT,*) ' -t, --tallies Write tally results from state point' write(OUTPUT_UNIT,*) ' -v, --version Show version information' write(OUTPUT_UNIT,*) ' -?, --help Show this message' @@ -1215,24 +1215,24 @@ contains if (entropy_on) then if (cmfd_run) then - message = " Bat./Gen. k(batch) Entropy Average k CMFD k CMFD Ent" + message = " Bat./Gen. k Entropy Average k CMFD k CMFD Ent" call write_message(1) message = " ========= ======== ======== ==================== ======== ========" call write_message(1) else - message = " Bat./Gen. k(batch) Entropy Average k" + message = " Bat./Gen. k Entropy Average k" call write_message(1) message = " ========= ======== ======== ====================" call write_message(1) end if else if (cmfd_run) then - message = " Bat./Gen. k(batch) Average k CMFD k" + message = " Bat./Gen. k Average k CMFD k" call write_message(1) message = " ========= ======== ==================== ========" call write_message(1) else - message = " Bat./Gen. k(batch) Average k" + message = " Bat./Gen. k Average k" call write_message(1) message = " ========= ======== ====================" call write_message(1) @@ -1251,11 +1251,17 @@ contains ! write out information about batch and generation write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(current_gen)) - write(UNIT=OUTPUT_UNIT, FMT='(11X)', ADVANCE='NO') + write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & + k_generation(overall_gen) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy(current_gen + gen_per_batch*(current_batch - 1)) + entropy(overall_gen) + + if (overall_gen - n_inactive*gen_per_batch > 1) then + write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & + keff, keff_std + end if ! next line write(UNIT=OUTPUT_UNIT, FMT=*) @@ -1273,14 +1279,14 @@ contains write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch)) write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_batch(current_batch) + k_generation(overall_gen) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & entropy(current_batch*gen_per_batch) ! write out accumulated k-effective if after first active batch - if (current_batch > n_inactive + 1) then + if (overall_gen - n_inactive*gen_per_batch > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & keff, keff_std else diff --git a/src/output_interface.F90 b/src/output_interface.F90 new file mode 100644 index 0000000000..5f421a6d25 --- /dev/null +++ b/src/output_interface.F90 @@ -0,0 +1,1127 @@ +module output_interface + + use constants + use error, only: warning + use global + use tally_header, only: TallyResult + +#ifdef HDF5 + use hdf5_interface +#elif MPI + use mpiio_interface +#endif + + implicit none + + ! Generic write procedure interface + interface write_data + module procedure write_double + module procedure write_double_1Darray + module procedure write_double_2Darray + module procedure write_double_3Darray + module procedure write_integer + module procedure write_integer_1Darray + module procedure write_integer_2Darray + module procedure write_integer_3Darray + module procedure write_long + module procedure write_string + end interface write_data + + ! Generic read procedure interface + interface read_data + module procedure read_double + module procedure read_double_1Darray + module procedure read_integer + module procedure read_integer_1Darray + module procedure read_long + module procedure read_string + end interface read_data + +contains + +!=============================================================================== +! FILE_CREATE creates a new file to write data to +!=============================================================================== + + subroutine file_create(filename, fh_str, proc_id) + + character(*), intent(in) :: filename ! name of file to be created + character(*), intent(in) :: fh_str ! parallel or serial HDF5 file + integer, optional, intent(in) :: proc_id ! processor rank to write from + + integer :: proc_create = 0 ! processor writing in serial (default master) + +#ifdef HDF5 +# ifdef MPI + ! Check for proc id + if (present(proc_id)) proc_create = proc_id + + ! Determine whether the file should be created by 1 or all procs + if (trim(fh_str) == 'serial') then + if(rank == proc_create) call hdf5_file_create(filename, hdf5_fh) + else + call hdf5_parallel_file_create(filename, hdf5_fh) + endif +# else + call hdf5_file_create(filename, hdf5_fh) +# endif +#elif MPI + call mpi_create_file(filename, mpi_fh) +#else + open(UNIT=UNIT_OUTPUT, FILE=filename, ACTION="write", & + STATUS='replace', ACCESS='stream') +#endif + + end subroutine file_create + +!=============================================================================== +! FILE_OPEN opens an existing file for reading or read/writing +!=============================================================================== + + subroutine file_open(filename, fh_str, mode, proc_id) + + character(*), intent(in) :: filename ! name of file to be opened + character(*), intent(in) :: fh_str ! parallel or serial HDF5 file + character(*), intent(in) :: mode ! file access mode + integer, optional, intent(in) :: proc_id ! processor rank to open file + + integer :: proc_open = 0 ! processor to open file (default is master) + +#ifdef HDF5 +# ifdef MPI + ! Check for proc_id + if (present(proc_id)) proc_open = proc_id + + ! Determine if the file should be opened by 1 or all procs + if (trim(fh_str) == 'serial') then + if (rank == proc_open) call hdf5_file_open(filename, hdf5_fh, mode) + else + call hdf5_parallel_file_open(filename, hdf5_fh, mode) + endif +# else + call hdf5_file_open(filename, hdf5_fh, mode) +# endif +#elif MPI + call mpi_open_file(filename, mpi_fh, mode) +#else + ! Check for read/write mode to open, default is read only + if (mode == 'w') then + open(UNIT=UNIT_OUTPUT, FILE=filename, ACTION='write', & + STATUS='old', ACCESS='stream') + else + open(UNIT=UNIT_OUTPUT, FILE=filename, ACTION='read', & + STATUS='old', ACCESS='stream') + end if +#endif + + end subroutine file_open + +!=============================================================================== +! FILE_CLOSE closes a file +!=============================================================================== + + subroutine file_close(fh_str, proc_id) + + character(*), intent(in) :: fh_str ! serial or parallel hdf5 file + integer, optional, intent(in) :: proc_id ! processor rank to close file + + integer :: proc_close = 0 ! processor to close file + +#ifdef HDF5 +# ifdef MPI + ! Check for proc_id + if (present(proc_id)) proc_close = proc_id + + ! Determine whether a file should be closed by 1 or all procs + if (trim(fh_str) == 'serial') then + if(rank == proc_close) call hdf5_file_close(hdf5_fh) + else + call hdf5_file_close(hdf5_fh) + endif +# else + call hdf5_file_close(hdf5_fh) +# endif +#elif MPI + call mpi_close_file(mpi_fh) +#else + close(UNIT=UNIT_OUTPUT) +#endif + + end subroutine file_close + +!=============================================================================== +! WRITE_DOUBLE writes double precision scalar data +!=============================================================================== + + subroutine write_double(buffer, name, group) + + real(8), intent(in) :: buffer ! data to write + character(*), intent(in) :: name ! name for data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_double(temp_group, name, buffer) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_write_double(mpi_fh, buffer) +#else + write(UNIT_OUTPUT) buffer +#endif + + end subroutine write_double + +!=============================================================================== +! READ_DOUBLE reads double precision scalar data +!=============================================================================== + + subroutine read_double(buffer, name, group, option) + + real(8), intent(inout) :: buffer ! read data to here + character(*), intent(in) :: name ! name for data + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in), optional :: option ! type of read + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif +# ifdef MPI + ! Check for option for reading default is independent + if (present(option)) then + if (option == 'collective') then + call hdf5_parallel_read_double(temp_group, name, buffer, & + H5FD_MPIO_COLLECTIVE_F) + else + call hdf5_parallel_read_double(temp_group, name, buffer, & + H5FD_MPIO_INDEPENDENT_F) + end if + else + ! Standard read call + call hdf5_read_double(temp_group, name, buffer) + end if +# else + ! Read the data serial + call hdf5_read_double(temp_group, name, buffer) +# endif + ! Check if HDf5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_read_double(mpi_fh, buffer) +#else + read(UNIT_OUTPUT) buffer +#endif + + end subroutine read_double + +!=============================================================================== +! WRITE_DOUBLE_1DARRAY writes double presicions 1-D array data +!=============================================================================== + + subroutine write_double_1Darray(buffer, name, group, length) + + integer, intent(in) :: length ! length of array to write + real(8), intent(in) :: buffer(:) ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_double_1Darray(temp_group, name, buffer, length) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_write_double_1Darray(mpi_fh, buffer, length) +#else + write(UNIT_OUTPUT) buffer(1:length) +#endif + + end subroutine write_double_1Darray + +!=============================================================================== +! READ_DOUBLE_1DARRAY reads double precision 1-D array data +!=============================================================================== + + subroutine read_double_1Darray(buffer, name, group, length, option) + + integer, intent(in) :: length ! length of array to read + real(8), intent(inout) :: buffer(:) ! read data to here + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in), optional :: option ! read option + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif +# ifdef MPI + ! Check for option for reading default is independent + if (present(option)) then + if (option == 'collective') then + call hdf5_parallel_read_double_1Darray(temp_group, name, buffer, & + length, H5FD_MPIO_COLLECTIVE_F) + else + call hdf5_parallel_read_double_1Darray(temp_group, name, buffer, & + length, H5FD_MPIO_INDEPENDENT_F) + end if + else + ! Standard read call + call hdf5_read_double_1Darray(temp_group, name, buffer, length) + end if +# else + ! Read the data serial + call hdf5_read_double_1Darray(temp_group, name, buffer, length) +# endif + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_read_double_1Darray(mpi_fh, buffer, length) +#else + read(UNIT_OUTPUT) buffer(1:length) +#endif + + end subroutine read_double_1Darray + +!=============================================================================== +! WRITE_DOUBLE_2DARRAY writes double precision 2-D array data +!=============================================================================== + + subroutine write_double_2Darray(buffer, name, group, length) + + integer, intent(in) :: length(2) ! dimension of array + real(8), intent(in) :: buffer(length(1),length(2)) ! the data + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_double_2Darray(temp_group, name, buffer, length) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#else + message = 'Double precision 2-D array writing not currently supported' + call warning() +#endif + + end subroutine write_double_2Darray + +!=============================================================================== +! WRITE_DOUBLE_3DARRAY writes double precision 3-D array data +!=============================================================================== + + subroutine write_double_3Darray(buffer, name, group, length) + + integer, intent(in) :: length(3) ! length of each dimension + real(8), intent(in) :: buffer(length(1),length(2),length(3)) + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_double_3Darray(temp_group, name, buffer, length) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#else + message = 'Double precision 3-D array writing not currently supported' + call warning() +#endif + + end subroutine write_double_3Darray + +!=============================================================================== +! WRITE_INTEGER writes integer scalar data +!=============================================================================== + + subroutine write_integer(buffer, name, group) + + integer, intent(in) :: buffer ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write data + call hdf5_write_integer(temp_group, name, buffer) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_write_integer(mpi_fh, buffer) +#else + write(UNIT_OUTPUT) buffer +#endif + + end subroutine write_integer + +!=============================================================================== +! READ_INTEGER reads integer scalar data +!=============================================================================== + + subroutine read_integer(buffer, name, group, option) + + integer, intent(inout) :: buffer ! read data to here + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in), optional :: option ! read option + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif +# ifdef MPI + ! Check for option for reading default is independent + if (present(option)) then + if (option == 'collective') then + call hdf5_parallel_read_integer(temp_group, name, buffer, & + H5FD_MPIO_COLLECTIVE_F) + else + call hdf5_parallel_read_integer(temp_group, name, buffer, & + H5FD_MPIO_INDEPENDENT_F) + end if + else + ! Standard read call + call hdf5_read_integer(temp_group, name, buffer) + end if +# else + ! Read the data serial + call hdf5_read_integer(temp_group, name, buffer) +# endif + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_read_integer(mpi_fh, buffer) +#else + read(UNIT_OUTPUT) buffer +#endif + + end subroutine read_integer + +!=============================================================================== +! WRITE_INTEGER_1DARRAY writes integer 1-D array data +!=============================================================================== + + subroutine write_integer_1Darray(buffer, name, group, length) + + integer, intent(in) :: length ! length of array to write + integer, intent(in) :: buffer(:) ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_integer_1Darray(temp_group, name, buffer, length) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_write_integer_1Darray(mpi_fh, buffer, length) +#else + write(UNIT_OUTPUT) buffer(1:length) +#endif + + end subroutine write_integer_1Darray + +!=============================================================================== +! READ_INTEGER_1DARRAY reads integer 1-D array data +!=============================================================================== + + subroutine read_integer_1Darray(buffer, name, group, length, option) + + integer, intent(in) :: length ! length of array to read + integer, intent(inout) :: buffer(:) ! read data to here + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in), optional :: option ! read option + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif +# ifdef MPI + ! Check for option for reading default is independent + if (present(option)) then + if (option == 'collective') then + call hdf5_parallel_read_integer_1Darray(temp_group, name, buffer, & + length, H5FD_MPIO_COLLECTIVE_F) + else + call hdf5_parallel_read_integer_1Darray(temp_group, name, buffer, & + length, H5FD_MPIO_INDEPENDENT_F) + end if + else + ! Standard read call + call hdf5_read_integer_1Darray(temp_group, name, buffer, length) + end if +# else + ! Read the data serial + call hdf5_read_integer_1Darray(temp_group, name, buffer, length) +# endif + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_read_integer_1Darray(mpi_fh, buffer, length) +#else + read(UNIT_OUTPUT) buffer(1:length) +#endif + + end subroutine read_integer_1Darray + +!=============================================================================== +! WRITE_INTEGER_2DARRAY writes integer 2-D array data +!=============================================================================== + + subroutine write_integer_2Darray(buffer, name, group, length) + + integer, intent(in) :: length(2) ! length of dimensions + integer, intent(in) :: buffer(length(1),length(2)) ! data + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_integer_2Darray(temp_group, name, buffer, length) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#else + message = 'Integer 2-D array writing not currently supported' + call warning() +#endif + + end subroutine write_integer_2Darray + +!=============================================================================== +! WRITE_INTEGER_3DARRAY writes integer 3-D array data +!=============================================================================== + + subroutine write_integer_3Darray(buffer, name, group, length) + + integer, intent(in) :: length(3) ! length of dimensions + integer, intent(in) :: buffer(length(1),length(2),length(3)) + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_integer_3Darray(temp_group, name, buffer, length) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#else + message = 'Integer 3-D array writing not currently supported' + call warning() +#endif + + end subroutine write_integer_3Darray + +!=============================================================================== +! WRITE_LONG writes long integer scalar data +!=============================================================================== + + subroutine write_long(buffer, name, group) + + integer(8), intent(in) :: buffer ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_long(temp_group, name, buffer, hdf5_integer8_t) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_write_long(mpi_fh, buffer) +#else + write(UNIT_OUTPUT) buffer +#endif + + end subroutine write_long + +!=============================================================================== +! READ_LONG reads long integer scalar data +!=============================================================================== + + subroutine read_long(buffer, name, group, option) + + integer(8), intent(inout) :: buffer ! read data to here + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in), optional :: option ! read option + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif +# ifdef MPI + ! Check for option for reading default is independent + if (present(option)) then + if (option == 'collective') then + call hdf5_parallel_read_long(temp_group, name, buffer, & + hdf5_integer8_t, H5FD_MPIO_COLLECTIVE_F) + else + call hdf5_parallel_read_long(temp_group, name, buffer, & + hdf5_integer8_t, H5FD_MPIO_INDEPENDENT_F) + end if + else + ! Standard read call + call hdf5_read_long(temp_group, name, buffer, hdf5_integer8_t) + end if +# else + ! Read the data serial + call hdf5_read_long(temp_group, name, buffer, hdf5_integer8_t) +# endif + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + call mpi_read_long(mpi_fh, buffer) +#else + read(UNIT_OUTPUT) buffer +#endif + + end subroutine read_long + +!=============================================================================== +! WRITE_STRING writes string data +!=============================================================================== + + subroutine write_string(buffer, name, group) + + character(*), intent(in) :: buffer ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + +#ifndef HDF5 +# ifdef MPI + integer :: n ! length of string buffer to write +# endif +#endif + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the data + call hdf5_write_string(temp_group, name, buffer, len(buffer)) + + ! Check if HDf5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + ! Length of string buffer to write + n = len(buffer) + + ! Write the data + call mpi_write_string(mpi_fh, buffer, n) +#else + write(UNIT_OUTPUT) buffer +#endif + + end subroutine write_string + +!=============================================================================== +! READ_STRING reads string data +!=============================================================================== + + subroutine read_string(buffer, name, group, option) + + character(*), intent(inout) :: buffer ! read data to here + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in), optional :: option ! read option + + integer :: n ! length of string to read to + + ! Length of string buffer to read + n = len(buffer) + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif +# ifdef MPI + ! Check for option for reading default is independent + if (present(option)) then + if (option == 'collective') then + call hdf5_parallel_read_string(temp_group, name, buffer, n, & + H5FD_MPIO_COLLECTIVE_F) + else + call hdf5_parallel_read_string(temp_group, name, buffer, n, & + H5FD_MPIO_INDEPENDENT_F) + end if + else + ! Standard read call + call hdf5_read_string(temp_group, name, buffer) + end if +# else + ! Read the data serial + call hdf5_read_string(temp_group, name, buffer) +# endif + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#elif MPI + + ! Read the data + call mpi_read_string(mpi_fh, buffer, n) +#else + read(UNIT_OUTPUT) buffer +#endif + + end subroutine read_string + +!=============================================================================== +! WRITE_ATTRIBUTE_STRING +!=============================================================================== + + subroutine write_attribute_string(var, attr_type, attr_str, group) + + character(*), intent(in) :: var ! variable name for attr + character(*), intent(in) :: attr_type ! attr identifier type + character(*), intent(in) :: attr_str ! string for attr id type + character(*), intent(in), optional :: group ! HDF5 group name + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + endif + + ! Write the attribute string + call hdf5_write_attribute_string(temp_group, var, attr_type, attr_str) + + ! Check if HDF5 group should be closed + if (present(group)) call hdf5_close_group() +#endif + + end subroutine write_attribute_string + +!=============================================================================== +! WRITE_TALLY_RESULT writes an OpenMC TallyResult type +!=============================================================================== + + subroutine write_tally_result(buffer, name, group, n1, n2) + + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in) :: name ! name of data + integer, intent(in) :: n1, n2 ! TallyResult dims + type(TallyResult), intent(in), target :: buffer(n1, n2) ! data to write + +#ifndef HDF5 +# ifndef MPI + integer :: j,k ! iteration counters +# endif +#endif + +#ifdef HDF5 + + ! Open up sub-group if present + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + end if + + ! Set overall size of vector to write + dims1(1) = n1*n2 + + ! Create up a dataspace for size + call h5screate_simple_f(1, dims1, dspace, hdf5_err) + + ! Create the dataset + call h5dcreate_f(temp_group, name, hdf5_tallyresult_t, dspace, dset, & + hdf5_err) + + ! Set pointer to first value and write + f_ptr = c_loc(buffer(1,1)) + call h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) + + ! Close ids + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + if (present(group)) then + call hdf5_close_group() + end if + +#elif MPI + + ! Write out tally buffer + call MPI_FILE_WRITE(mpi_fh, buffer, n1*n2, MPI_TALLYRESULT, & + MPI_STATUS_IGNORE, mpiio_err) + +#else + + ! Write out tally buffer + do k = 1, n2 + do j = 1, n1 + write(UNIT_OUTPUT) buffer(j,k) % sum + write(UNIT_OUTPUT) buffer(j,k) % sum_sq + end do + end do + +#endif + + end subroutine write_tally_result + +!=============================================================================== +! READ_TALLY_RESULT reads OpenMC TallyResult data +!=============================================================================== + + subroutine read_tally_result(buffer, name, group, n1, n2) + + character(*), intent(in), optional :: group ! HDF5 group name + character(*), intent(in) :: name ! name of data + integer, intent(in) :: n1, n2 ! TallyResult dims + type(TallyResult), intent(inout), target :: buffer(n1, n2) ! read data here + +#ifndef HDF5 +# ifndef MPI + integer :: j,k ! iteration counters +# endif +#endif + +#ifdef HDF5 + + ! Open up sub-group if present + if (present(group)) then + call hdf5_open_group(group) + else + temp_group = hdf5_fh + end if + + ! Open the dataset + call h5dopen_f(temp_group, name, dset, hdf5_err) + + ! Set pointer to first value and write + f_ptr = c_loc(buffer(1,1)) + call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) + + ! Close ids + call h5dclose_f(dset, hdf5_err) + if (present(group)) call hdf5_close_group() + +#elif MPI + + ! Write out tally buffer + call MPI_FILE_READ(mpi_fh, buffer, n1*n2, MPI_TALLYRESULT, & + MPI_STATUS_IGNORE, mpiio_err) + +#else + + ! Write out tally buffer + do k = 1, n2 + do j = 1, n1 + read(UNIT_OUTPUT) buffer(j,k) % sum + read(UNIT_OUTPUT) buffer(j,k) % sum_sq + end do + end do + +#endif + + end subroutine read_tally_result + + +!=============================================================================== +! WRITE_SOURCE_BANK writes OpenMC source_bank data +!=============================================================================== + + subroutine write_source_bank() + +#ifdef HDF5 + integer(8) :: offset(1) ! source data offset +#elif MPI + integer(MPI_OFFSET_KIND) :: offset ! offset of data + integer :: size_offset_kind ! the data offset kind + integer :: size_bank ! size of bank to write +#endif + +#ifdef HDF5 +# ifdef MPI + + ! Set size of total dataspace for all procs and rank + dims1(1) = n_particles + hdf5_rank = 1 + + ! Create that dataspace + call h5screate_simple_f(hdf5_rank, dims1, dspace, hdf5_err) + + ! Create the dataset for that dataspace + call h5dcreate_f(hdf5_fh, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err) + + ! Close the dataspace + call h5sclose_f(dspace, hdf5_err) + + ! Create another data space but for each proc individually + dims1(1) = work + call h5screate_simple_f(hdf5_rank, dims1, memspace, hdf5_err) + + ! Get the individual local proc dataspace + call h5dget_space_f(dset, dspace, hdf5_err) + + ! Select hyperslab for this dataspace + offset(1) = bank_first - 1_8 + call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims1, 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(1)) + + ! 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 + + ! Set size + dims1(1) = work + hdf5_rank = 1 + + ! Create dataspace + call h5screate_simple_f(hdf5_rank, dims1, dspace, hdf5_err) + + ! Create dataset + call h5dcreate_f(hdf5_fh, "source_bank", hdf5_bank_t, & + dspace, dset, hdf5_err) + + ! Set up pointer to data + f_ptr = c_loc(source_bank(1)) + + ! Write dataset to file + call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err) + + ! Close all ids + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + +# endif + +#elif MPI + + ! Get current offset for master + if (master) call MPI_FILE_GET_POSITION(mpi_fh, offset, mpiio_err) + + ! Determine offset on master process and broadcast to all processors + call MPI_SIZEOF(offset, size_offset_kind, mpi_err) + select case (size_offset_kind) + case (4) + call MPI_BCAST(offset, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpi_err) + case (8) + call MPI_BCAST(offset, 1, MPI_INTEGER8, 0, MPI_COMM_WORLD, mpi_err) + end select + + ! Set the proper offset for source data on this processor + call MPI_TYPE_SIZE(MPI_BANK, size_bank, mpi_err) + offset = offset + size_bank*maxwork*rank + + ! Write all source sites + call MPI_FILE_WRITE_AT(mpi_fh, offset, source_bank(1), work, MPI_BANK, & + MPI_STATUS_IGNORE, mpiio_err) + +#else + + ! Write out source sites + write(UNIT_OUTPUT) source_bank + +#endif + + end subroutine write_source_bank + +!=============================================================================== +! READ_SOURCE_BANK reads OpenMC source_bank data +!=============================================================================== + + subroutine read_source_bank() + +#ifdef HDF5 + integer(8) :: offset(1) ! offset of data +#elif MPI + integer(MPI_OFFSET_KIND) :: offset ! offset of data + integer :: size_offset_kind ! the data offset kind + integer :: size_bank ! size of bank to read +#endif + +#ifdef HDF5 +# ifdef MPI + + ! Set size of total dataspace for all procs and rank + dims1(1) = n_particles + hdf5_rank = 1 + + ! Open the dataset + call h5dopen_f(hdf5_fh, "source_bank", dset, hdf5_err) + + ! Create another data space but for each proc individually + dims1(1) = work + call h5screate_simple_f(hdf5_rank, dims1, memspace, hdf5_err) + + ! Get the individual local proc dataspace + call h5dget_space_f(dset, dspace, hdf5_err) + + ! Select hyperslab for this dataspace + offset(1) = bank_first - 1_8 + call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims1, 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(1)) + + ! Read data from file in parallel + call h5dread_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 + + ! Open dataset + call h5dopen_f(hdf5_fh, "source_bank", dset, hdf5_err) + + ! Set up pointer to data + f_ptr = c_loc(source_bank(1)) + + ! Read dataset from file + call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err) + + ! Close all ids + call h5dclose_f(dset, hdf5_err) + +# endif + +#elif MPI + + ! Get current offset for master + if (master) call MPI_FILE_GET_POSITION(mpi_fh, offset, mpiio_err) + + ! Determine offset on master process and broadcast to all processors + call MPI_SIZEOF(offset, size_offset_kind, mpi_err) + select case (size_offset_kind) + case (4) + call MPI_BCAST(offset, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpi_err) + case (8) + call MPI_BCAST(offset, 1, MPI_INTEGER8, 0, MPI_COMM_WORLD, mpi_err) + end select + + ! Set the proper offset for source data on this processor + call MPI_TYPE_SIZE(MPI_BANK, size_bank, mpi_err) + offset = offset + size_bank*maxwork*rank + + ! Write all source sites + call MPI_FILE_READ_AT(mpi_fh, offset, source_bank(1), work, MPI_BANK, & + MPI_STATUS_IGNORE, mpiio_err) + +#else + + ! Write out source sites + read(UNIT_OUTPUT) source_bank + +#endif + + end subroutine read_source_bank + +end module output_interface diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 07c7f5475e..7c732800fa 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -7,7 +7,7 @@ module particle_restart use geometry_header, only: BASE_UNIVERSE use global use particle_header, only: deallocate_coord - use output, only: write_message + use output, only: write_message, print_particle use physics, only: transport use random_lcg, only: set_particle_seed use source, only: initialize_particle @@ -38,8 +38,6 @@ contains subroutine read_hdf5_particle_restart() - integer(HSIZE_T) :: dims1(1) - ! write meessage message = "Loading particle restart file " // trim(path_particle_restart) & // "..." @@ -53,8 +51,8 @@ contains call hdf5_read_integer(hdf5_particle_file, 'current_batch', current_batch) call hdf5_read_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch) call hdf5_read_integer(hdf5_particle_file, 'current_gen', current_gen) - call hdf5_read_long(hdf5_particle_file, 'n_particles', n_particles) - call hdf5_read_long(hdf5_particle_file, 'id', p % id) + call hdf5_read_long(hdf5_particle_file, 'n_particles', n_particles, hdf5_integer8_t) + call hdf5_read_long(hdf5_particle_file, 'id', p % id, hdf5_integer8_t) call hdf5_read_double(hdf5_particle_file, 'weight', p % wgt) call hdf5_read_double(hdf5_particle_file, 'energy', p % E) dims1 = (/3/) @@ -81,6 +79,9 @@ contains subroutine read_binary_particle_restart() + integer :: filetype + integer :: revision + ! write meessage message = "Loading particle restart file " // trim(path_particle_restart) & // "..." @@ -91,6 +92,8 @@ contains ACCESS='stream') ! read data from file + read(UNIT_PARTICLE) filetype + read(UNIT_PARTICLE) revision read(UNIT_PARTICLE) current_batch read(UNIT_PARTICLE) gen_per_batch read(UNIT_PARTICLE) current_gen @@ -142,11 +145,7 @@ contains call transport() ! write output if particle made it - write(ou,*) 'Particle Successfully Transport:' - write(ou,*) 'WEIGHT:', p % wgt - write(ou,*) 'ENERGY:', p % E - write(ou,*) 'LOCATION:', p % coord % xyz - write(ou,*) 'ANGLE:', p % coord % uvw + call print_particle() end subroutine run_particle_restart diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 index 8b2e7a9012..12cd17c34b 100644 --- a/src/particle_restart_write.F90 +++ b/src/particle_restart_write.F90 @@ -47,7 +47,6 @@ contains subroutine write_particle_restart_hdf5() character(MAX_FILE_LEN) :: filename - integer(HSIZE_T) :: dims1(1) type(Bank), pointer :: src => null() ! set up file name @@ -61,11 +60,15 @@ contains src => source_bank(current_work) ! write data to file + call hdf5_write_integer(hdf5_particle_file, 'filetype', & + FILETYPE_PARTICLE_RESTART) + call hdf5_write_integer(hdf5_particle_file, 'revision', & + REVISION_PARTICLE_RESTART) call hdf5_write_integer(hdf5_particle_file, 'current_batch', current_batch) call hdf5_write_integer(hdf5_particle_file, 'gen_per_batch', gen_per_batch) call hdf5_write_integer(hdf5_particle_file, 'current_gen', current_gen) - call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles) - call hdf5_write_long(hdf5_particle_file, 'id', p % id) + call hdf5_write_long(hdf5_particle_file, 'n_particles', n_particles, hdf5_integer8_t) + call hdf5_write_long(hdf5_particle_file, 'id', p % id, hdf5_integer8_t) call hdf5_write_double(hdf5_particle_file, 'weight', src % wgt) call hdf5_write_double(hdf5_particle_file, 'energy', src % E) dims1 = (/3/) @@ -102,6 +105,8 @@ contains src => source_bank(current_work) ! write data to file + write(UNIT_PARTICLE) FILETYPE_PARTICLE_RESTART + write(UNIT_PARTICLE) REVISION_PARTICLE_RESTART write(UNIT_PARTICLE) current_batch write(UNIT_PARTICLE) gen_per_batch write(UNIT_PARTICLE) current_gen diff --git a/src/source.F90 b/src/source.F90 index 4075422366..9ffea50089 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -9,7 +9,6 @@ module source use particle_header, only: deallocate_coord use physics, only: maxwell_spectrum, watt_spectrum use random_lcg, only: prn, set_particle_seed - use state_point, only: read_source_binary use string, only: to_str #ifdef MPI @@ -38,7 +37,8 @@ contains ! Read the source from a binary file instead of sampling from some ! assumed source distribution - call read_source_binary() + message = 'This feature is currently disabled and will be added back in.' + call fatal_error() else ! Generation source sites from specified distribution in user input @@ -167,8 +167,7 @@ contains p % id = bank_first + index_source - 1 ! set random number seed - particle_seed = ((current_batch - 1)*gen_per_batch + & - current_gen - 1)*n_particles + p % id + particle_seed = (overall_gen - 1)*n_particles + p % id call set_particle_seed(particle_seed) ! set particle trace diff --git a/src/state_point.F90 b/src/state_point.F90 index cf11676f28..4628221581 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -1,564 +1,375 @@ module state_point !=============================================================================== -! STATE_POINT -- This module handles writing and reading binary state point +! STATE_POINT -- This module handles writing and reading state point ! files. State points are contain complete tally results, source sites, and ! various other data. They can be used to restart a run or to reconstruct ! confidence intervals for tallies (this requires post-processing via Python ! scripts). ! -! Modifications to this module should be made with care. There are essentially -! three different ways to write or read state points: 1) normal Fortran file -! I/O, 2) MPI-IO, and 3) HDF5. The HDF5 functionality is contained in the -! hdf5_interface module. If you plan to change the state point, you will need to -! change all methods. You should also increment REVISION_STATEPOINT in the -! constants module. -! ! State points can be written at any batch during a simulation, or at specified ! intervals, using the tag. !=============================================================================== - use error, only: warning, fatal_error - use global - use math, only: t_percentile - use output, only: write_message, print_batch_keff, time_stamp - use string, only: to_str - use tally_header, only: TallyObject -#ifdef MPI - use mpi -#endif + use constants + use error, only: fatal_error, warning + use global + use output, only: write_message, time_stamp + use string, only: to_str + use output_interface + use tally_header, only: TallyObject implicit none contains !=============================================================================== -! WRITE_STATE_POINT creates a state point binary file that can be used for -! restarting a run or for getting intermediate tally results +! WRITE_STATE_POINT !=============================================================================== subroutine write_state_point() - integer :: i ! loop index -#ifdef MPI - integer :: fh ! file handle - integer :: n ! temporary array length - integer :: temp ! temporary variable - integer :: size_offset_kind ! size of MPI_OFFSET_KIND (bytes) - integer :: size_bank ! size of MPI_BANK type - integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) -#else - integer :: j, k ! loop indices -#endif character(MAX_FILE_LEN) :: filename + integer :: i + integer :: j + integer, allocatable :: temp_array(:) type(TallyObject), pointer :: t => null() - ! Set filename for binary state point + ! Set filename for state point filename = trim(path_output) // 'statepoint.' // & - trim(to_str(current_batch)) // '.binary' + trim(to_str(current_batch)) + + ! Append appropriate extension +#ifdef HDF5 + filename = trim(filename) // '.h5' +#else + filename = trim(filename) // '.binary' +#endif ! Write message message = "Creating state point " // trim(filename) // "..." call write_message(1) -#ifdef MPI - ! ========================================================================== - ! PARALLEL I/O USING MPI-2 ROUTINES + ! Create statepoint file + call file_create(filename, 'serial') - ! Open binary source file for reading - call MPI_FILE_OPEN(MPI_COMM_WORLD, filename, MPI_MODE_CREATE + & - MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err) + if (master) then + ! Write file type + call write_data(FILETYPE_STATEPOINT, "filetype") - ! ========================================================================== - ! RUN INFORMATION AND TALLY METADATA + ! Write revision number for state point file + call write_data(REVISION_STATEPOINT, "revision") - if (master) call write_state_point_header(fh) + ! Write OpenMC version + call write_data(VERSION_MAJOR, "version_major") + call write_data(VERSION_MINOR, "version_minor") + call write_data(VERSION_RELEASE, "version_release") - ! ========================================================================== - ! TALLY RESULTS + ! Write current date and time + call write_data(time_stamp(), "date_and_time") + ! Write path to input + call write_data(path_input, "path") + + ! Write out random number seed + call write_data(seed, "seed") + + ! Write run information + call write_data(run_mode, "run_mode") + call write_data(n_particles, "n_particles") + call write_data(n_batches, "n_batches") + + ! Write out current batch number + call write_data(current_batch, "current_batch") + + ! Write out information for eigenvalue run + if (run_mode == MODE_EIGENVALUE) then + call write_data(n_inactive, "n_inactive") + call write_data(gen_per_batch, "gen_per_batch") + call write_data(k_generation, "k_generation", & + length=current_batch*gen_per_batch) + call write_data(entropy, "entropy", length=current_batch*gen_per_batch) + call write_data(k_col_abs, "k_col_abs") + call write_data(k_col_tra, "k_col_tra") + call write_data(k_abs_tra, "k_abs_tra") + call write_data(k_combined, "k_combined", length=2) + end if + + ! Write number of meshes + call write_data(n_meshes, "n_meshes", group="tallies") + + ! Write information for meshes + MESH_LOOP: do i = 1, n_meshes + call write_data(meshes(i) % id, "id", & + group="tallies/mesh" // to_str(i)) + call write_data(meshes(i) % type, "type", & + group="tallies/mesh" // to_str(i)) + call write_data(meshes(i) % n_dimension, "n_dimension", & + group="tallies/mesh" // to_str(i)) + call write_data(meshes(i) % dimension, "dimension", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension) + call write_data(meshes(i) % lower_left, "lower_left", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension) + call write_data(meshes(i) % upper_right, "upper_right", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension) + call write_data(meshes(i) % width, "width", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension) + end do MESH_LOOP + + ! Write number of tallies + call write_data(n_tallies, "n_tallies", group="tallies") + + ! Write all tally information except results + TALLY_METADATA: do i = 1, n_tallies + !Get pointer to tally + t => tallies(i) + + ! Write id + call write_data(t % id, "id", group="tallies/tally" // to_str(i)) + + ! Write number of realizations + call write_data(t % n_realizations, "n_realizations", & + group="tallies/tally" // to_str(i)) + + ! Write size of each tally + call write_data(t % total_score_bins, "total_score_bins", & + group="tallies/tally" // to_str(i)) + call write_data(t % total_filter_bins, "total_filter_bins", & + group="tallies/tally" // to_str(i)) + + ! Write number of filters + call write_data(t % n_filters, "n_filters", & + group="tallies/tally" // to_str(i)) + + ! Write filter information + FILTER_LOOP: do j = 1, t % n_filters + + ! Write type of filter + call write_data(t % filters(j) % type, "type", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j)) + + ! Write number of bins for this filter + call write_data(t % filters(j) % n_bins, "n_bins", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j)) + + ! Write bins + if (t % filters(j) % type == FILTER_ENERGYIN .or. & + t % filters(j) % type == FILTER_ENERGYOUT) then + call write_data(t % filters(j) % real_bins, "bins", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), & + length=size(t % filters(j) % real_bins)) + else + call write_data(t % filters(j) % int_bins, "bins", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), & + length=size(t % filters(j) % int_bins)) + end if + + end do FILTER_LOOP + + ! Write number of nuclide bins + call write_data(t % n_nuclide_bins, "n_nuclide_bins", & + group="tallies/tally" // to_str(i)) + + ! Set up nuclide bin array and then write + allocate(temp_array(t % n_nuclide_bins)) + NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins + if (t % nuclide_bins(j) > 0) then + temp_array(j) = nuclides(t % nuclide_bins(j)) % zaid + else + temp_array(j) = t % nuclide_bins(j) + end if + end do NUCLIDE_LOOP + call write_data(temp_array, "nuclide_bins", & + group="tallies/tally" // to_str(i), length=t % n_nuclide_bins) + deallocate(temp_array) + + ! Write number of score bins, score bins, and scatt order + call write_data(t % n_score_bins, "n_score_bins", & + group="tallies/tally" // to_str(i)) + call write_data(t % score_bins, "score_bins", & + group="tallies/tally" // to_str(i), length=t % n_score_bins) + call write_data(t % scatt_order, "scatt_order", & + group="tallies/tally" // to_str(i), length=t % n_score_bins) + + ! Write number of user score bins + call write_data(t % n_user_score_bins, "n_user_score_bins", & + group="tallies/tally" // to_str(i)) + + end do TALLY_METADATA + + end if + + ! Check for the no-tally-reduction method if (.not. reduce_tallies) then ! If using the no-tally-reduction method, we need to collect tally ! results before writing them to the state point file. - call write_tally_results_nr(fh) + call write_tally_results_nr() elseif (master) then - ! Write number of realizations - call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + + ! Write number of global realizations + call write_data(n_realizations, "n_realizations") ! Write global tallies - call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, global_tallies, N_GLOBAL_TALLIES, & - MPI_TALLYRESULT, MPI_STATUS_IGNORE, mpi_err) + call write_data(N_GLOBAL_TALLIES, "n_global_tallies") + call write_tally_result(global_tallies, "global_tallies", & + n1=N_GLOBAL_TALLIES, n2=1) + ! Write tallies if (tallies_on) then + ! Indicate that tallies are on - temp = 1 - call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call write_data(1, "tallies_present", group="tallies") ! Write all tally results TALLY_RESULTS: do i = 1, n_tallies + + ! Set point to current tally t => tallies(i) - n = size(t % results, 1) * size(t % results, 2) - call MPI_FILE_WRITE(fh, t % results, n, MPI_TALLYRESULT, & - MPI_STATUS_IGNORE, mpi_err) + ! Write sum and sum_sq for each bin + call write_tally_result(t % results, "results", & + group="tallies/tally" // to_str(i), & + n1=size(t % results, 1), n2=size(t % results, 2)) + end do TALLY_RESULTS + else - ! Indicate that tallies are off - temp = 0 - call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + + ! Indicate tallies are off + call write_data(0, "tallies_present", group="tallies") + end if + end if - ! ========================================================================== - ! SOURCE BANK + ! Check for eigenvalue calculation + if (run_mode == MODE_EIGENVALUE .and. source_write) then - if (run_mode == MODE_EIGENVALUE) then + ! Check for writing source out separately if (source_separate) then - ! If the user has specified that the source sites should be written in - ! a separate file, we make a call to the appropriate subroutine to - ! write it separately - if (source_write) then - path_source = "source." // trim(to_str(current_batch)) // ".binary" - call write_source_binary() - end if - elseif (source_write) then - ! Otherwise, write the source sites in the state point file - - ! Get current offset for master - if (master) call MPI_FILE_GET_POSITION(fh, offset, mpi_err) - - ! Determine offset on master process and broadcast to all processors - call MPI_SIZEOF(offset, size_offset_kind, mpi_err) - select case (size_offset_kind) - case (4) - call MPI_BCAST(offset, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpi_err) - case (8) - call MPI_BCAST(offset, 1, MPI_INTEGER8, 0, MPI_COMM_WORLD, mpi_err) - end select - - ! Set proper offset for source data on this processor - call MPI_TYPE_SIZE(MPI_BANK, size_bank, mpi_err) - offset = offset + size_bank*maxwork*rank - - ! Write all source sites - call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, & - MPI_STATUS_IGNORE, mpi_err) - end if - end if - - ! Close binary source file - call MPI_FILE_CLOSE(fh, mpi_err) + ! Close statepoint file + call file_close('serial') + ! Set filename for source + filename = trim(path_output) // 'source.' // & + trim(to_str(current_batch)) +#ifdef HDF5 + filename = trim(filename) // '.h5' #else - ! Open binary state point file for writing - open(UNIT=UNIT_STATE, FILE=filename, STATUS='replace', & - ACCESS='stream') - - ! Write revision number for state point file - write(UNIT_STATE) REVISION_STATEPOINT - - ! Write OpenMC version - write(UNIT_STATE) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE - - ! Write current date and time - write(UNIT_STATE) time_stamp() - - ! Write path to input - write(UNIT_STATE) path_input - - ! Write out random number seed - write(UNIT_STATE) seed - - ! Write run information - write(UNIT_STATE) run_mode, n_particles, n_batches - - ! Write out current batch number - write(UNIT_STATE) current_batch - - ! Write out information for eigenvalue run - if (run_mode == MODE_EIGENVALUE) then - write(UNIT_STATE) n_inactive, gen_per_batch - write(UNIT_STATE) k_batch(1:current_batch) - write(UNIT_STATE) entropy(1:current_batch*gen_per_batch) - write(UNIT_STATE) k_col_abs - write(UNIT_STATE) k_col_tra - write(UNIT_STATE) k_abs_tra - write(UNIT_STATE) k_combined - end if - - ! Write number of meshes - write(UNIT_STATE) n_meshes - - ! Write information for meshes - MESH_LOOP: do i = 1, n_meshes - write(UNIT_STATE) meshes(i) % id - write(UNIT_STATE) meshes(i) % type - write(UNIT_STATE) meshes(i) % n_dimension - write(UNIT_STATE) meshes(i) % dimension - write(UNIT_STATE) meshes(i) % lower_left - write(UNIT_STATE) meshes(i) % upper_right - write(UNIT_STATE) meshes(i) % width - end do MESH_LOOP - - ! Write number of tallies - write(UNIT_STATE) n_tallies - - TALLY_METADATA: do i = 1, n_tallies - ! Get pointer to tally - t => tallies(i) - - ! Write id - write(UNIT_STATE) t % id - - ! Number of realizations - write(UNIT_STATE) t % n_realizations - - ! Write size of each dimension of tally results array - write(UNIT_STATE) t % total_score_bins - write(UNIT_STATE) t % total_filter_bins - - ! Write number of filters - write(UNIT_STATE) t % n_filters - - FILTER_LOOP: do j = 1, t % n_filters - ! Write type of filter - write(UNIT_STATE) t % filters(j) % type - - ! Write number of bins for this filter - write(UNIT_STATE) t % filters(j) % n_bins - - ! Write filter bins - if (t % filters(j) % type == FILTER_ENERGYIN .or. & - t % filters(j) % type == FILTER_ENERGYOUT) then - write(UNIT_STATE) t % filters(j) % real_bins - else - write(UNIT_STATE) t % filters(j) % int_bins - end if - end do FILTER_LOOP - - ! Write number of nuclide bins - write(UNIT_STATE) t % n_nuclide_bins - - ! Write nuclide bins - NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins - if (t % nuclide_bins(j) > 0) then - write(UNIT_STATE) nuclides(t % nuclide_bins(j)) % zaid - else - write(UNIT_STATE) t % nuclide_bins(j) - end if - end do NUCLIDE_LOOP - - ! Write number of score bins, score bins, and scatt order - write(UNIT_STATE) t % n_score_bins - write(UNIT_STATE) t % score_bins - write(UNIT_STATE) t % scatt_order - - ! Write number of user score bins - write(UNIT_STATE) t % n_user_score_bins - end do TALLY_METADATA - - ! Number of realizations for global tallies - write(UNIT_STATE) n_realizations - - ! Write out global tallies sum and sum_sq - write(UNIT_STATE) N_GLOBAL_TALLIES - GLOBAL_TALLIES_LOOP: do i = 1, N_GLOBAL_TALLIES - write(UNIT_STATE) global_tallies(i) % sum - write(UNIT_STATE) global_tallies(i) % sum_sq - end do GLOBAL_TALLIES_LOOP - - if (tallies_on) then - ! Indicate that tallies are on - write(UNIT_STATE) 1 - - TALLY_RESULTS: do i = 1, n_tallies - ! Get pointer to tally - t => tallies(i) - - ! Write tally sum and sum_sq for each bin - do k = 1, size(t % results, 2) - do j = 1, size(t % results, 1) - write(UNIT_STATE) t % results(j,k) % sum - write(UNIT_STATE) t % results(j,k) % sum_sq - end do - end do - end do TALLY_RESULTS - else - ! Indicate that tallies are off - write(UNIT_STATE) 0 - end if - - ! Write out source bank - if (run_mode == MODE_EIGENVALUE) then - if (source_separate) then - ! If the user has specified that the source sites should be written in - ! a separate file, we make a call to the appropriate subroutine to - ! write it separately - - if (source_write) then - path_source = "source." // trim(to_str(current_batch)) // ".binary" - call write_source_binary() - end if - elseif (source_write) then - ! Otherwise, write the source sites in the state point file - - write(UNIT_STATE) source_bank - end if - end if - - ! Close binary state point file - close(UNIT_STATE) + filename = trim(filename) // '.binary' #endif + ! Write message + message = "Creating source file " // trim(filename) // "..." + call write_message(1) + + ! Create statepoint file + call file_create(filename, 'parallel') + +#ifdef HDF5 +# ifdef MPI + else + ! Close HDF5 serial file and reopen in parallel + call file_close('serial') + call file_open(filename, 'parallel', 'w') +# endif +#endif + + end if + + ! Write out source + call write_source_bank() + + ! Close file, all files in parallel mode + call file_close('parallel') ! even if no MPI, this will work for HDF5 + + else + + ! Close file if not in eigenvalue mode or no source writing + call file_close('serial') + + end if + end subroutine write_state_point -#ifdef MPI -!=============================================================================== -! WRITE_STATE_POINT_HEADER uses MPI-IO routines to write basic run information -! and tally metadata -!=============================================================================== - - subroutine write_state_point_header(fh) - - integer, intent(inout) :: fh ! file handle - - integer :: i ! loop index - integer :: j ! loop index - integer :: n ! temporary array length - type(TallyObject), pointer :: t => null() - - ! Write revision number for state point file - call MPI_FILE_WRITE(fh, REVISION_STATEPOINT, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write OpenMC version - call MPI_FILE_WRITE(fh, VERSION_MAJOR, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, VERSION_MINOR, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, VERSION_RELEASE, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write current date and time - call MPI_FILE_WRITE(fh, time_stamp(), 19, MPI_CHARACTER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write path to input - call MPI_FILE_WRITE(fh, path_input, MAX_FILE_LEN, MPI_CHARACTER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write out random number seed - call MPI_FILE_WRITE(fh, seed, 1, MPI_INTEGER8, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write run information - call MPI_FILE_WRITE(fh, run_mode, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, n_particles, 1, MPI_INTEGER8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, n_batches, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write out current batch number - call MPI_FILE_WRITE(fh, current_batch, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write out information for eigenvalue run - if (run_mode == MODE_EIGENVALUE) then - call MPI_FILE_WRITE(fh, n_inactive, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, gen_per_batch, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, k_batch, current_batch, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, entropy, current_batch*gen_per_batch, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, k_col_abs, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, k_col_tra, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, k_abs_tra, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, k_combined, 2, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - end if - - ! Write number of meshes - call MPI_FILE_WRITE(fh, n_meshes, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write information for meshes - MESH_LOOP: do i = 1, n_meshes - call MPI_FILE_WRITE(fh, meshes(i) % id, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, meshes(i) % type, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, meshes(i) % n_dimension, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - n = meshes(i) % n_dimension - call MPI_FILE_WRITE(fh, meshes(i) % dimension, n, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, meshes(i) % lower_left, n, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, meshes(i) % upper_right, n, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, meshes(i) % width, n, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - end do MESH_LOOP - - ! Write number of tallies - call MPI_FILE_WRITE(fh, n_tallies, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - TALLY_METADATA: do i = 1, n_tallies - ! Get pointer to tally - t => tallies(i) - - ! Write id - call MPI_FILE_WRITE(fh, t % id, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write number of realizations - call MPI_FILE_WRITE(fh, t % n_realizations, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write size of each tally - call MPI_FILE_WRITE(fh, t % total_score_bins, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, t % total_filter_bins, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write number of filters - call MPI_FILE_WRITE(fh, t % n_filters, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - FILTER_LOOP: do j = 1, t % n_filters - ! Write type of filter - call MPI_FILE_WRITE(fh, t % filters(j) % type, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write number of bins for this filter - call MPI_FILE_WRITE(fh, t % filters(j) % n_bins, & - 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - - ! Write bins - if (t % filters(j) % type == FILTER_ENERGYIN .or. & - t % filters(j) % type == FILTER_ENERGYOUT) then - n = size(t % filters(j) % real_bins) - call MPI_FILE_WRITE(fh, t % filters(j) % real_bins, n, & - MPI_REAL8, MPI_STATUS_IGNORE, mpi_err) - else - n = size(t % filters(j) % int_bins) - call MPI_FILE_WRITE(fh, t % filters(j) % int_bins, n, & - MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - end if - end do FILTER_LOOP - - ! Write number of nuclide bins - call MPI_FILE_WRITE(fh, t % n_nuclide_bins, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - - ! Write nuclide bins - NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins - if (t % nuclide_bins(j) > 0) then - call MPI_FILE_WRITE(fh, nuclides(t % nuclide_bins(j)) % zaid, & - 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - else - call MPI_FILE_WRITE(fh, t % nuclide_bins(j), 1, & - MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - end if - end do NUCLIDE_LOOP - - ! Write number of score bins, score bins, and scatt order - call MPI_FILE_WRITE(fh, t % n_score_bins, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, t % score_bins, t % n_score_bins, & - MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_WRITE(fh, t % scatt_order, t % n_score_bins, & - MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - - ! Write number of user score bins - call MPI_FILE_WRITE(fh, t % n_user_score_bins, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - end do TALLY_METADATA - - end subroutine write_state_point_header -#endif - -#ifdef MPI !=============================================================================== ! WRITE_TALLY_RESULTS_NR !=============================================================================== - subroutine write_tally_results_nr(fh) - - integer, intent(in) :: fh ! file handle + subroutine write_tally_results_nr() integer :: i ! loop index integer :: n ! number of filter bins integer :: m ! number of score bins - integer :: temp ! temporary variable integer :: n_bins ! total number of bins real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results - real(8) :: global_temp(2,N_GLOBAL_TALLIES) + real(8), target :: global_temp(2,N_GLOBAL_TALLIES) real(8) :: dummy ! temporary receive buffer for non-root reduces type(TallyObject), pointer :: t => null() + type(TallyResult), allocatable :: tallyresult_temp(:,:) ! ========================================================================== ! COLLECT AND WRITE GLOBAL TALLIES if (master) then ! Write number of realizations - call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call write_data(n_realizations, "n_realizations") ! Write number of global tallies - call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call write_data(N_GLOBAL_TALLIES, "n_global_tallies") end if ! Copy global tallies into temporary array for reducing n_bins = 2 * N_GLOBAL_TALLIES global_temp(1,:) = global_tallies(:) % sum - global_temp(2,:) = global_tallies(:) % sum_sq + global_temp(2,:) = global_tallies(:) % sum_sq if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into a ! receive buffer without having a temporary variable +#ifdef MPI call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, & 0, MPI_COMM_WORLD, mpi_err) - - ! Write out global tallies sum and sum_sq - call MPI_FILE_WRITE(fh, global_temp, n_bins, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) +#endif ! Transfer values to value on master if (current_batch == n_batches) then global_tallies(:) % sum = global_temp(1,:) global_tallies(:) % sum_sq = global_temp(2,:) end if + + ! Put reduced value in temporary tally result + allocate(tallyresult_temp(N_GLOBAL_TALLIES, 1)) + tallyresult_temp(:,1) % sum = global_temp(1,:) + tallyresult_temp(:,1) % sum_sq = global_temp(2,:) + + + ! Write out global tallies sum and sum_sq + call write_tally_result(tallyresult_temp, "global_tallies", & + n1=N_GLOBAL_TALLIES, n2=1) + + ! Deallocate temporary tally result + deallocate(tallyresult_temp) else ! Receive buffer not significant at other processors +#ifdef MPI call MPI_REDUCE(global_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & 0, MPI_COMM_WORLD, mpi_err) +#endif end if if (tallies_on) then ! Indicate that tallies are on if (master) then - temp = 1 - call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call write_data(1, "tallies_present", group="tallies") end if ! Write all tally results @@ -579,23 +390,34 @@ contains if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable +#ifdef MPI call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) - - ! Write reduced tally results to file - call MPI_FILE_WRITE(fh, tally_temp, n_bins, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - +#endif ! At the end of the simulation, store the results back in the ! regular TallyResults array if (current_batch == n_batches) then t % results(:,:) % sum = tally_temp(1,:,:) t % results(:,:) % sum_sq = tally_temp(2,:,:) end if + + ! Put in temporary tally result + allocate(tallyresult_temp(m,n)) + tallyresult_temp(:,:) % sum = tally_temp(1,:,:) + tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:) + + ! Write reduced tally results to file + call write_tally_result(t % results, "results", & + group="tallies/tally" // to_str(i), n1=m, n2=n) + + ! Deallocate temporary tally result + deallocate(tallyresult_temp) else ! Receive buffer not significant at other processors +#ifdef MPI call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & 0, MPI_COMM_WORLD, mpi_err) +#endif end if ! Deallocate temporary copy of tally results @@ -604,684 +426,300 @@ contains else if (master) then ! Indicate that tallies are off - temp = 0 - call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call write_data(0, "tallies_present", group="tallies") end if end if end subroutine write_tally_results_nr -#endif !=============================================================================== -! LOAD_STATE_POINT loads data from a state point file to either continue a run -! or to print intermediate tally results +! LOAD_STATE_POINT !=============================================================================== subroutine load_state_point() - integer :: i, j ! loop indices - integer :: mode ! specified run mode - integer :: temp(3) ! temporary variable - integer, allocatable :: int_array(:) - real(8), allocatable :: real_array(:) - character(19) :: current_time ! current date and time + character(MAX_FILE_LEN) :: filename character(MAX_FILE_LEN) :: path_temp - -#ifdef MPI - integer :: fh ! file handle - integer :: n ! temporary array size - integer :: size_offset_kind ! size of MPI_OFFSET_KIND (bytes) - integer :: size_bank ! size of MPI_BANK type - integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) -#else - integer :: k ! loop index -#endif + character(19) :: current_time + integer :: i + integer :: j + integer :: int_array(3) + integer, allocatable :: temp_array(:) + real(8) :: real_array(3) + type(TallyObject), pointer :: t => null() ! Write message message = "Loading state point " // trim(path_state_point) // "..." call write_message(1) -#ifdef MPI - ! Open binary state point file for reading - call MPI_FILE_OPEN(MPI_COMM_WORLD, path_state_point, MPI_MODE_RDONLY, & - MPI_INFO_NULL, fh, mpi_err) + ! Open file for reading + call file_open(path_state_point, 'parallel', 'r') - ! ========================================================================== - ! RUN INFORMATION AND TALLY METADATA + ! Read filetype + call read_data(int_array(1), "filetype", option="collective") - ! Raad revision number for state point file and make sure it matches with + ! Read revision number for state point file and make sure it matches with ! current version - call MPI_FILE_READ_ALL(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= REVISION_STATEPOINT) then - message = "State point binary version does not match current version " & - // "in OpenMC." + call read_data(int_array(1), "revision", option="collective") + if (int_array(1) /= REVISION_STATEPOINT) then + message = "State point version does not match current version " & + // "in OpenMC." call fatal_error() end if ! Read OpenMC version - call MPI_FILE_READ_ALL(fh, temp, 3, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= VERSION_MAJOR .or. temp(2) /= VERSION_MINOR & - .or. temp(3) /= VERSION_RELEASE) then - message = "State point file was created with a different version " // & - "of OpenMC." + call read_data(int_array(1), "version_major", option="collective") + call read_data(int_array(2), "version_minor", option="collective") + call read_data(int_array(3), "version_release", option="collective") + if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR & + .or. int_array(3) /= VERSION_RELEASE) then + message = "State point file was created with a different version " & + // "of OpenMC." call warning() end if ! Read date and time - call MPI_FILE_READ_ALL(fh, current_time, 19, MPI_CHARACTER, & - MPI_STATUS_IGNORE, mpi_err) + call read_data(current_time, "date_and_time", option="collective") ! Read path to input - call MPI_FILE_READ_ALL(fh, path_temp, MAX_FILE_LEN, MPI_CHARACTER, & - MPI_STATUS_IGNORE, mpi_err) + call read_data(path_temp, "path", option="collective") ! Read and overwrite random number seed - call MPI_FILE_READ_ALL(fh, seed, 1, MPI_INTEGER8, & - MPI_STATUS_IGNORE, mpi_err) + call read_data(seed, "seed", option="collective") - ! Read and overwrite run information - call MPI_FILE_READ_ALL(fh, mode, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, n_particles, 1, MPI_INTEGER8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, n_batches, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + ! Read and overwrite run information except number of batches + call read_data(run_mode, "run_mode", option="collective") + call read_data(n_particles, "n_particles", option="collective") + call read_data(int_array(1), "n_batches", option="collective") + + ! Take maximum of statepoint n_batches and input n_batches + n_batches = max(n_batches, int_array(1)) ! Read batch number to restart at - call MPI_FILE_READ_ALL(fh, restart_batch, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call read_data(restart_batch, "current_batch", option="collective") ! Read information specific to eigenvalue run - if (mode == MODE_EIGENVALUE) then - call MPI_FILE_READ_ALL(fh, n_inactive, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, gen_per_batch, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, k_batch, restart_batch, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, entropy, restart_batch*gen_per_batch, & - MPI_REAL8, MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, k_col_abs, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, k_col_tra, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ_ALL(fh, k_abs_tra, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - allocate(real_array(2)) - call MPI_FILE_READ_ALL(fh, real_array, 2, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - deallocate(real_array) + if (run_mode == MODE_EIGENVALUE) then + call read_data(int_array(1), "n_inactive", option="collective") + call read_data(gen_per_batch, "gen_per_batch", option="collective") + call read_data(k_generation, "k_generation", & + length=restart_batch*gen_per_batch, option="collective") + call read_data(entropy, "entropy", length=restart_batch*gen_per_batch, & + option="collective") + call read_data(k_col_abs, "k_col_abs", option="collective") + call read_data(k_col_tra, "k_col_tra", option="collective") + call read_data(k_abs_tra, "k_abs_tra", option="collective") + call read_data(real_array(1:2), "k_combined", length=2, & + option="collective") + + ! Take maximum of statepoint n_inactive and input n_inactive + n_inactive = max(n_inactive, int_array(1)) end if - if (master) then - ! Read number of meshes - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= n_meshes) then - message = "Number of meshes does not match in state point." + ! Read number of meshes + call read_data(n_meshes, "n_meshes", group="tallies", option="collective") + + ! Read and overwrite mesh information + MESH_LOOP: do i = 1, n_meshes + call read_data(meshes(i) % id, "id", & + group="tallies/mesh" // to_str(i), option="collective") + call read_data(meshes(i) % type, "type", & + group="tallies/mesh" // to_str(i), option="collective") + call read_data(meshes(i) % n_dimension, "n_dimension", & + group="tallies/mesh" // to_str(i), option="collective") + call read_data(meshes(i) % dimension, "dimension", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension, option="collective") + call read_data(meshes(i) % lower_left, "lower_left", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension, option="collective") + call read_data(meshes(i) % upper_right, "upper_right", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension, option="collective") + call read_data(meshes(i) % width, "width", & + group="tallies/mesh" // to_str(i), & + length=meshes(i) % n_dimension, option="collective") + end do MESH_LOOP + + ! Read and overwrite number of tallies + call read_data(n_tallies, "n_tallies", group="tallies", option="collective") + + ! Read in tally metadata + TALLY_METADATA: do i = 1, n_tallies + + ! Get pointer to tally + t => tallies(i) + + ! Read tally id + call read_data(t % id, "id", group="tallies/tally" // to_str(i), & + option="collective") + + ! Read number of realizations + call read_data(t % n_realizations, "n_realizations", & + group="tallies/tally" // to_str(i), option="collective") + + ! Read size of tally results + call read_data(int_array(1), "total_score_bins", & + group="tallies/tally" // to_str(i), option="collective") + call read_data(int_array(2), "total_filter_bins", & + group="tallies/tally" // to_str(i), option="collective") + + ! Check size of tally results array + if (int_array(1) /= t % total_score_bins .and. & + int_array(2) /= t % total_filter_bins) then + message = "Input file tally structure is different from restart." call fatal_error() end if - MESH_LOOP: do i = 1, n_meshes - ! Read id, mesh type, and dimension - call MPI_FILE_READ(fh, temp, 3, MPI_INTEGER, MPI_STATUS_IGNORE, & - mpi_err) + ! Read number of filters + call read_data(t % n_filters, "n_filters", & + group="tallies/tally" // to_str(i), option="collective") - ! Skip mesh data - call MPI_FILE_GET_POSITION(fh, offset, mpi_err) - offset = offset + temp(3)*(4 + 3*8) - call MPI_FILE_SEEK(fh, offset, MPI_SEEK_SET, mpi_err) - end do MESH_LOOP + ! Read filter information + FILTER_LOOP: do j = 1, t % n_filters - ! Read number of tallies and make sure it matches - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= n_tallies) then - message = "Number of tallies does not match in state point." - call fatal_error() - end if + ! Read type of filter + call read_data(t % filters(j) % type, "type", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), & + option="collective") - TALLY_METADATA: do i = 1, n_tallies - ! Read tally id - call MPI_FILE_READ(fh, tallies(i) % id, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + ! Read number of bins for this filter + call read_data(t % filters(j) % n_bins, "n_bins", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), & + option="collective") - ! Read number of realizations for global tallies - call MPI_FILE_READ(fh, tallies(i) % n_realizations, 1, & - MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - - ! Read dimensions of tally filters and results and make sure they - ! match - call MPI_FILE_READ(fh, temp, 2, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= size(tallies(i) % results, 1) .or. & - temp(2) /= size(tallies(i) % results, 2)) then - message = "Tally dimensions do not match in state point." - call fatal_error() + ! Read bins + if (t % filters(j) % type == FILTER_ENERGYIN .or. & + t % filters(j) % type == FILTER_ENERGYOUT) then + call read_data(t % filters(j) % real_bins, "bins", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), & + length=size(t % filters(j) % real_bins), option="collective") + else + call read_data(t % filters(j) % int_bins, "bins", & + group="tallies/tally" // trim(to_str(i)) // "/filter" // to_str(j), & + length=size(t % filters(j) % int_bins), option="collective") end if - ! Read number of filters - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + end do FILTER_LOOP - FILTER_LOOP: do j = 1, temp(1) - ! Read filter type and number of bins - call MPI_FILE_READ(fh, temp(2), 2, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + ! Read number of nuclide bins + call read_data(t % n_nuclide_bins, "n_nuclide_bins", & + group="tallies/tally" // to_str(i), option="collective") - ! Read filter bins - select case (temp(2)) - case (FILTER_MESH) - allocate(int_array(1)) - call MPI_FILE_READ(fh, int_array, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - deallocate(int_array) - case (FILTER_ENERGYIN, FILTER_ENERGYOUT) - allocate(real_array(temp(3) + 1)) - call MPI_FILE_READ(fh, real_array, temp(3) + 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpi_err) - deallocate(real_array) - case default - allocate(int_array(temp(3))) - call MPI_FILE_READ(fh, int_array, temp(3), MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - deallocate(int_array) - end select - end do FILTER_LOOP + ! Set up nuclide bin array and then write + allocate(temp_array(t % n_nuclide_bins)) + call read_data(temp_array, "nuclide_bins", & + group="tallies/tally" // to_str(i), length=t % n_nuclide_bins, & + option="collective") + NUCLIDE_LOOP: do j = 1, t % n_nuclide_bins + if (temp_array(j) > 0) then + nuclides(t % nuclide_bins(j)) % zaid = temp_array(j) + else + t % nuclide_bins(j) = temp_array(j) + end if + end do NUCLIDE_LOOP + deallocate(temp_array) - ! Read number of nuclides - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + ! Write number of score bins, score bins, and scatt order + call read_data(t % n_score_bins, "n_score_bins", & + group="tallies/tally" // to_str(i), option="collective") + call read_data(t % score_bins, "score_bins", & + group="tallies/tally" // to_str(i), length=t % n_score_bins, & + option="collective") + call read_data(t % scatt_order, "scatt_order", & + group="tallies/tally" // to_str(i), length=t % n_score_bins, & + option="collective") - ! Read nuclide bins - allocate(int_array(temp(1))) - call MPI_FILE_READ(fh, int_array, temp(1), MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - deallocate(int_array) + ! Write number of user score bins + call read_data(t % n_user_score_bins, "n_user_score_bins", & + group="tallies/tally" // to_str(i), option="collective") - ! Read number of score bins, score bins, and scatt_order - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - allocate(int_array(temp(1))) - call MPI_FILE_READ(fh, int_array, temp(1), MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - call MPI_FILE_READ(fh, int_array, temp(1), MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - deallocate(int_array) - - ! Read number of user score bins - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) - end do TALLY_METADATA + end do TALLY_METADATA + + ! Read tallies to master + if (master) then ! Read number of realizations for global tallies - call MPI_FILE_READ(fh, n_realizations, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpi_err) + call read_data(n_realizations, "n_realizations") - ! Read number of global tallies and make sure it matches - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= N_GLOBAL_TALLIES) then + ! Read number of global tallies + call read_data(int_array(1), "n_global_tallies") + if (int_array(1) /= N_GLOBAL_TALLIES) then message = "Number of global tallies does not match in state point." call fatal_error() end if ! Read global tally data - call MPI_FILE_READ(fh, global_tallies, N_GLOBAL_TALLIES, & - MPI_TALLYRESULT, MPI_STATUS_IGNORE, mpi_err) + call read_tally_result(global_tallies, "global_tallies", & + n1=N_GLOBAL_TALLIES, n2=1) ! Check if tally results are present - call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) + call read_data(int_array(1), "tallies_present", group="tallies") - ! ======================================================================= - ! TALLY RESULTS - - ! Read sum and sum squared - if (temp(1) == 1) then + ! Read in sum and sum squared + if (int_array(1) == 1) then TALLY_RESULTS: do i = 1, n_tallies - n = size(tallies(i) % results, 1) * size(tallies(i) % results, 2) - call MPI_FILE_READ(fh, tallies(i) % results, n, MPI_TALLYRESULT, & - MPI_STATUS_IGNORE, mpi_err) + + ! Set pointer to tally + t => tallies(i) + + ! Read sum and sum_sq for each bin + call read_tally_result(t % results, "results", & + group="tallies/tally" // to_str(i), & + n1=size(t % results, 1), n2=size(t % results, 2)) + end do TALLY_RESULTS end if end if - ! ========================================================================== - ! SOURCE BANK + ! Read source if in eigenvalue mode + if (run_mode == MODE_EIGENVALUE .and. run_mode /= MODE_TALLIES) then - if (run_mode == MODE_EIGENVALUE) then - ! Get current offset for master - if (master) call MPI_FILE_GET_POSITION(fh, offset, mpi_err) + ! Check if source was written out separately + if (source_separate) then - ! Determine offset on master process and broadcast to all processors - call MPI_SIZEOF(offset, size_offset_kind, mpi_err) - select case (size_offset_kind) - case (4) - call MPI_BCAST(offset, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpi_err) - case (8) - call MPI_BCAST(offset, 1, MPI_INTEGER8, 0, MPI_COMM_WORLD, mpi_err) - end select - - ! Set proper offset for source data on this processor - call MPI_TYPE_SIZE(MPI_BANK, size_bank, mpi_err) - offset = offset + size_bank*maxwork*rank - - ! Write all source sites - call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, & - MPI_STATUS_IGNORE, mpi_err) - end if - - ! Close binary state point file - call MPI_FILE_CLOSE(fh, mpi_err) + ! Close statepoint file + call file_close('parallel') + ! Set filename for source + filename = trim(path_output) // 'source.' // & + trim(to_str(restart_batch)) +#ifdef HDF5 + filename = trim(filename) // '.h5' #else - ! Open binary state point file for writing - open(UNIT=UNIT_STATE, FILE=path_state_point, STATUS='old', & - ACCESS='stream') - - ! Raad revision number for state point file and make sure it matches with - ! current version - read(UNIT_STATE) temp(1) - if (temp(1) /= REVISION_STATEPOINT) then - message = "State point binary version does not match current version " & - // "in OpenMC." - call fatal_error() - end if - - ! Read OpenMC version - read(UNIT_STATE) temp(1:3) - if (temp(1) /= VERSION_MAJOR .or. temp(2) /= VERSION_MINOR & - .or. temp(3) /= VERSION_RELEASE) then - message = "State point file was created with a different version " // & - "of OpenMC." - call warning() - end if - - ! Read date and time - read(UNIT_STATE) current_time - - ! Read path - read(UNIT_STATE) path_temp - - ! Read and overwrite random number seed - read(UNIT_STATE) seed - - ! Read and overwrite run information - read(UNIT_STATE) mode, n_particles, n_batches - - ! Read batch number to restart at - read(UNIT_STATE) restart_batch - - ! Read information specific to eigenvalue run - if (mode == MODE_EIGENVALUE) then - read(UNIT_STATE) n_inactive, gen_per_batch - read(UNIT_STATE) k_batch(1:restart_batch) - read(UNIT_STATE) entropy(1:restart_batch*gen_per_batch) - read(UNIT_STATE) k_col_abs - read(UNIT_STATE) k_col_tra - read(UNIT_STATE) k_abs_tra - allocate(real_array(2)) - read(UNIT_STATE) real_array - deallocate(real_array) - end if - - if (master) then - ! Read number of meshes - read(UNIT_STATE) temp(1) - if (temp(1) /= n_meshes) then - message = "Number of meshes does not match in state point." - call fatal_error() - end if - - MESH_LOOP: do i = 1, n_meshes - ! Read id, mesh type, and dimension - read(UNIT_STATE) temp(1:3) - - ! Allocate temporary arrays - allocate(int_array(temp(3))) - allocate(real_array(temp(3))) - - ! Read dimension, lower_left, upper_right, width - read(UNIT_STATE) int_array - read(UNIT_STATE) real_array - read(UNIT_STATE) real_array - read(UNIT_STATE) real_array - - ! Deallocate temporary arrays - deallocate(int_array) - deallocate(real_array) - end do MESH_LOOP - - ! Read number of tallies and make sure it matches - read(UNIT_STATE) temp(1) - if (temp(1) /= n_tallies) then - message = "Number of tallies does not match in state point." - call fatal_error() - end if - - TALLY_METADATA: do i = 1, n_tallies - ! Read id - read(UNIT_STATE) temp(1) - - ! Read number of realizations - read(UNIT_STATE) tallies(i) % n_realizations - - ! Read dimensions of tally filters and results and make sure they - ! match - read(UNIT_STATE) temp(1:2) - if (temp(1) /= size(tallies(i) % results, 1) .or. & - temp(2) /= size(tallies(i) % results, 2)) then - message = "Tally dimensions do not match in state point." - call fatal_error() - end if - - ! Read number of filters - read(UNIT_STATE) temp(1) - - FILTER_LOOP: do j = 1, temp(1) - ! Read filter type and number of bins - read(UNIT_STATE) temp(2:3) - - ! Read filter bins - select case (temp(2)) - case (FILTER_MESH) - allocate(int_array(1)) - read(UNIT_STATE) int_array - deallocate(int_array) - case (FILTER_ENERGYIN, FILTER_ENERGYOUT) - allocate(real_array(temp(3) + 1)) - read(UNIT_STATE) real_array - deallocate(real_array) - case default - allocate(int_array(temp(3))) - read(UNIT_STATE) int_array - deallocate(int_array) - end select - end do FILTER_LOOP - - ! Read number of nuclides - read(UNIT_STATE) temp(1) - - ! Read nuclide bins - allocate(int_array(temp(1))) - read(UNIT_STATE) int_array - deallocate(int_array) - - ! Read number of results - read(UNIT_STATE) temp(1) - - ! Read results bins and scatt_order - allocate(int_array(temp(1))) - read(UNIT_STATE) int_array - read(UNIT_STATE) int_array - deallocate(int_array) - - ! Read number of user bins - read(UNIT_STATE) temp(1) - end do TALLY_METADATA - - ! Read number of realizations for global tallies - read(UNIT_STATE) n_realizations - - ! Read number of global tallies and make sure it matches - read(UNIT_STATE) temp(1) - if (temp(1) /= N_GLOBAL_TALLIES) then - message = "Number of global tallies does not match in state point." - call fatal_error() - end if - - ! Read global tally data - do i = 1, N_GLOBAL_TALLIES - read(UNIT_STATE) global_tallies(i) % sum - read(UNIT_STATE) global_tallies(i) % sum_sq - end do - - ! Read sum and sum squared - read(UNIT_STATE) temp(1) - if (temp(1) == 1) then - TALLY_RESULTS: do i = 1, n_tallies - do k = 1, size(tallies(i) % results, 2) - do j = 1, size(tallies(i) % results, 1) - read(UNIT_STATE) tallies(i) % results(j,k) % sum - read(UNIT_STATE) tallies(i) % results(j,k) % sum_sq - end do - end do - end do TALLY_RESULTS - end if - end if - - ! Read source bank for eigenvalue run - if (mode == MODE_EIGENVALUE .and. run_mode /= MODE_TALLIES) then - read(UNIT_STATE) source_bank - end if - - ! Close binary state point file - close(UNIT_STATE) + filename = trim(filename) // '.binary' #endif + ! Write message + message = "Loading source file " // trim(filename) // "..." + call write_message(1) + + ! Create statepoint file + call file_open(filename, 'parallel', 'r') + + end if + + ! Write out source + call read_source_bank() + + ! Close file + if (source_separate) then + call file_close('parallel') + else + call file_close('parallel') + end if + + else + + ! Close file if not in eigenvalue mode + call file_close('parallel') + + end if + end subroutine load_state_point -!=============================================================================== -! REPLAY_BATCH_HISTORY displays batch keff and entropy for each batch stored in -! a state point file -!=============================================================================== - - subroutine replay_batch_history - - integer :: n = 0 ! number of realizations - real(8), save :: temp(2) = ZERO ! temporary values for keff - real(8) :: alpha ! significance level for CI - real(8) :: t_value ! t-value for confidence intervals - - ! Write message at beginning - if (current_batch == 1) then - message = "Replaying history from state point..." - call write_message(1) - end if - - ! Add to number of realizations - if (current_batch > n_inactive) then - n = n + 1 - - temp(1) = temp(1) + k_batch(current_batch) - temp(2) = temp(2) + k_batch(current_batch)*k_batch(current_batch) - - ! calculate mean keff - keff = temp(1) / n - - if (n > 1) then - if (confidence_intervals) then - ! Calculate t-value for confidence intervals - alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n - 1) - else - t_value = ONE - end if - - keff_std = t_value * sqrt((temp(2)/n - keff*keff)/(n - 1)) - end if - else - keff = k_batch(current_batch) - end if - - ! print out batch keff - if (master) call print_batch_keff() - - ! Write message at end - if (current_batch == restart_batch) then - message = "Resuming simulation..." - call write_message(1) - end if - - end subroutine replay_batch_history - -!=============================================================================== -! WRITE_SOURCE writes out the final source distribution to a binary file that -! can be used as a starting source in a new simulation -!=============================================================================== - - subroutine write_source_binary() - -#ifdef MPI - integer :: fh ! file handle - integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) - - ! ========================================================================== - ! PARALLEL I/O USING MPI-2 ROUTINES - - ! Open binary source file for reading - call MPI_FILE_OPEN(MPI_COMM_WORLD, path_source, MPI_MODE_CREATE + & - MPI_MODE_WRONLY, MPI_INFO_NULL, fh, mpi_err) - - if (master) then - offset = 0 - call MPI_FILE_WRITE_AT(fh, offset, n_particles, 1, MPI_INTEGER8, & - MPI_STATUS_IGNORE, mpi_err) - end if - - ! Set proper offset for source data on this processor - offset = 8*(1 + rank*maxwork*8) - - ! Write all source sites - call MPI_FILE_WRITE_AT(fh, offset, source_bank(1), work, MPI_BANK, & - MPI_STATUS_IGNORE, mpi_err) - - ! Close binary source file - call MPI_FILE_CLOSE(fh, mpi_err) - -#else - ! ========================================================================== - ! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES - - ! Open binary source file for writing - open(UNIT=UNIT_SOURCE, FILE=path_source, STATUS='replace', & - ACCESS='stream') - - ! Write the number of particles - write(UNIT=UNIT_SOURCE) n_particles - - ! Write information from the source bank - write(UNIT=UNIT_SOURCE) source_bank(1:work) - - ! Close binary source file - close(UNIT=UNIT_SOURCE) -#endif - - end subroutine write_source_binary - -!=============================================================================== -! READ_SOURCE_BINARY reads a source distribution from a source.binary file and -! initializes the source bank -!=============================================================================== - - subroutine read_source_binary() - - integer :: i ! loop over repeating sites - integer(8) :: n_sites ! number of sites in binary file - integer :: n_repeat ! number of times to repeat a site -#ifdef MPI - integer :: fh ! file handle - integer(MPI_OFFSET_KIND) :: offset ! offset in memory (0=beginning of file) - integer :: n_read ! number of sites to read on a single process -#endif - -#ifdef MPI - ! ========================================================================== - ! PARALLEL I/O USING MPI-2 ROUTINES - - ! Open binary source file for reading - call MPI_FILE_OPEN(MPI_COMM_WORLD, path_source, MPI_MODE_RDONLY, & - MPI_INFO_NULL, fh, mpi_err) - - ! Read number of source sites in file - offset = 0 - call MPI_FILE_READ_AT(fh, offset, n_sites, 1, MPI_INTEGER8, & - MPI_STATUS_IGNORE, mpi_err) - - if (n_particles > n_sites) then - ! Determine number of sites to read and offset - if (rank <= mod(n_sites,int(n_procs,8)) - 1) then - n_read = int(n_sites/n_procs) + 1 - offset = 8*(1 + rank*n_read*8) - else - n_read = int(n_sites/n_procs) - offset = 8*(1 + (rank*n_read + mod(n_sites,int(n_procs,8)))*8) - end if - - ! Read source sites - call MPI_FILE_READ_AT(fh, offset, source_bank(1), n_read, MPI_BANK, & - MPI_STATUS_IGNORE, mpi_err) - - ! Let's say we have 30 sites and we need to fill in 200. This do loop - ! will fill in sites 31 - 180. - - n_repeat = int(work / n_read) - do i = 1, n_repeat - 1 - source_bank(i*n_read + 1:(i+1)*n_read) = & - source_bank((i-1)*n_read + 1:i*n_read) - end do - - ! This final statement would fill sites 181 - 200 in the above example. - - if (mod(work, int(n_repeat*n_read,8)) > 0) then - source_bank(n_repeat*n_read + 1:work) = & - source_bank(1:work - n_repeat * n_read) - end if - - else - ! Set proper offset for source data on this processor - offset = 8*(1 + rank*maxwork*8) - - ! Read all source sites - call MPI_FILE_READ_AT(fh, offset, source_bank(1), work, MPI_BANK, & - MPI_STATUS_IGNORE, mpi_err) - - ! Close binary source file - call MPI_FILE_CLOSE(fh, mpi_err) - end if - -#else - ! ========================================================================== - ! SERIAL I/O USING FORTRAN INTRINSIC ROUTINES - - ! Open binary source file for reading - open(UNIT=UNIT_SOURCE, FILE=path_source, STATUS='old', & - ACCESS='stream') - - ! Read number of source sites in file - read(UNIT=UNIT_SOURCE) n_sites - - if (n_particles > n_sites) then - ! The size of the source file is smaller than the number of particles we - ! need. Thus, read all sites and then duplicate sites as necessary. - - read(UNIT=UNIT_SOURCE) source_bank(1:n_sites) - - ! Let's say we have 300 sites and we need to fill in 1000. This do loop - ! will fill in sites 301 - 900. - - n_repeat = int(n_particles / n_sites) - do i = 1, n_repeat - 1 - source_bank(i*n_sites + 1:(i+1)*n_sites) = & - source_bank((i-1)*n_sites + 1:i*n_sites) - end do - - ! This final statement would fill sites 901 - 1000 in the above example. - - source_bank(n_repeat*n_sites + 1:n_particles) = & - source_bank(1:n_particles - n_repeat * n_sites) - - else - ! The size of the source file is bigger than or equal to the number of - ! particles we need for one generation. Thus, we can just read as many - ! sites as we need. - - read(UNIT=UNIT_SOURCE) source_bank(1:n_particles) - - end if - - ! Close binary source file - close(UNIT=UNIT_SOURCE) -#endif - - end subroutine read_source_binary - + subroutine read_source +! TODO write this routine +! TODO what if n_particles does not match source bank + end subroutine read_source end module state_point diff --git a/src/tally.F90 b/src/tally.F90 index 3aa3e34855..b7b55189ad 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1757,12 +1757,6 @@ contains end do if (run_mode == MODE_EIGENVALUE) then - ! Get the current batch estimate of k_analog for displaying to output - ! --- this has to be performed after reduce_tally_values and before - ! accumulate_result - - k_batch(current_batch) = global_tallies(K_TRACKLENGTH) % value - if (active_batches) then ! Accumulate products of different estimators of k k_col = global_tallies(K_COLLISION) % value / total_weight diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 50b2297876..917ba19e18 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -147,8 +147,11 @@ class StatePoint(object): self._read_metadata() def _read_metadata(self): + # Read filetype + self.filetype = self._get_int(path='filetype')[0] + # Read statepoint revision - self.revision = self._get_int(path='revision_statepoint')[0] + self.revision = self._get_int(path='revision')[0] # Read OpenMC version if self._hdf5: @@ -179,7 +182,8 @@ class StatePoint(object): if self.run_mode == 2: self.n_inactive = self._get_int(path='n_inactive')[0] self.gen_per_batch = self._get_int(path='gen_per_batch')[0] - self.k_batch = self._get_double(self.current_batch, path='k_batch') + self.k_batch = self._get_double( + self.current_batch*self.gen_per_batch, path='k_generation') self.entropy = self._get_double( self.current_batch*self.gen_per_batch, path='entropy') if self.revision >= 8: