From 52233d0173494368ebe87210e748cf5403a382c3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Jul 2013 15:10:50 -0400 Subject: [PATCH] Implimented HDF5 support. --- docs/source/usersguide/processing.rst | 2 +- src/DEPENDENCIES | 1 + src/Makefile | 4 +- src/particle_track.F90 | 136 +++++++++++++++----------- src/utils/track.py | 43 +++++--- 5 files changed, 112 insertions(+), 74 deletions(-) diff --git a/docs/source/usersguide/processing.rst b/docs/source/usersguide/processing.rst index bd0d19a1a4..c12041214a 100644 --- a/docs/source/usersguide/processing.rst +++ b/docs/source/usersguide/processing.rst @@ -373,7 +373,7 @@ The settings.xml file can dictate that specific particle tracks are output. The 1 2 4 -After running OpenMC, the directory should contain a file of the form "track_(batch #)_(generation #)_(particle #).binary" for each particle tracked. These binary track files can be converted into VTK poly data files with the "track.py" utility. The usage of track.py is of the form "track.py [-o OUT] IN" where OUT is the optional output filename and IN is one or more filenames describing track binary files. The default output name is "track.pvtp". A common usage of track.py is "track.py track*.binary" which will use the data from all binary track files in the directory to write a "track.pvtp" VTK output file. The .pvtp file can then be read and plotted by 3d visualization programs such as Paraview. +After running OpenMC, the directory should contain a file of the form "track_(batch #)_(generation #)_(particle #).(binary or h5)" for each particle tracked. These track files can be converted into VTK poly data files with the "track.py" utility. The usage of track.py is of the form "track.py [-o OUT] IN" where OUT is the optional output filename and IN is one or more filenames describing track files. The default output name is "track.pvtp". A common usage of track.py is "track.py track*.binary" which will use the data from all binary track files in the directory to write a "track.pvtp" VTK output file. The .pvtp file can then be read and plotted by 3d visualization programs such as Paraview. diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index ea81f9ce9f..064acd707c 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -275,6 +275,7 @@ particle_header.o: constants.o particle_track.o: constants.o particle_track.o: global.o +particle_track.o: output_interface.o particle_track.o: string.o particle_restart.o: bank_header.o diff --git a/src/Makefile b/src/Makefile index c4a9c3dbef..05cb7df6e8 100644 --- a/src/Makefile +++ b/src/Makefile @@ -21,7 +21,7 @@ DEBUG = no PROFILE = no OPTIMIZE = no MPI = no -HDF5 = yes +HDF5 = no PETSC = no #=============================================================================== @@ -29,7 +29,7 @@ PETSC = no #=============================================================================== MPI_DIR = /opt/mpich/3.0.4-$(COMPILER) -HDF5_DIR = /usr/local/hdf5 +HDF5_DIR = /opt/hdf5/1.8.11-$(COMPILER) PHDF5_DIR = /opt/phdf5/1.8.11-$(COMPILER) PETSC_DIR = /opt/petsc/3.3-p6-$(COMPILER) diff --git a/src/particle_track.F90 b/src/particle_track.F90 index bc11497fbd..852a07e100 100644 --- a/src/particle_track.F90 +++ b/src/particle_track.F90 @@ -1,110 +1,136 @@ +!=============================================================================== +! PARTICLE_TRACK handles output of particle tracks to disk. +! +! TODO: This module writes binary output files via a stream access i.e. it +! writes the particle's location to disk everytime write_particle_track is +! called. But HDF5 does not allow for a simple implimentation of stream access. +! An attempt was made to write HDF5 files in a stream fashion using an +! extendable dataset, but it always failed at the write call because of an out- +! of-bounds error. The extendable dataset code is commented out by enclosing it +! within "#if 0"/"#endif" blocks, and has been replaced by code that writes +! particle coordinates to an ever-expanding array which is written to disk in +! its entirety when finalize_particle_track is called. +! +! In the interest of consistency, maybe this module should be changed so that it +! uses the extendable HDF5 dataset, or it writes binary output all in one call +! using the output_interface (like the current HDF5 implimentation). +!=============================================================================== + module particle_track use constants use global - !use output_interface only: file_create, file_close, write_double1Darray - use string, only: to_str + use string, only: to_str -#ifdef HDF5 +#if 0 use hdf5 #endif +#ifdef HDF5 + use output_interface, only: file_create, file_close, write_data +#endif implicit none +#if 0 + integer, private :: hdf5_err ! HDF error code + integer(HID_T), private :: track_dset ! dataset handle + integer(HID_T), private :: track_dspace ! dataspace handle + integer(HID_T), private :: track_fh ! HDF file handle + integer(HID_T), private :: cparms ! chunk parameters + integer(HSIZE_T), private :: n_tracks ! total number of tracks +#endif #ifdef HDF5 - integer, private :: hdf5_err ! HDF5 error code - integer(HID_T), private :: track_dset ! dataset handle - integer(HID_T), private :: track_dspace ! master dataspace handle - integer(HID_T), private :: track_selspace ! selected dataspace handle - integer(HID_T), private :: track_fh ! HDF5 file handle - integer(HID_T), private :: cparms ! chunk parameters - integer(HSIZE_T), private :: n_tracks ! total number of tracks + character(MAX_FILE_LEN), private :: fname ! file name + integer, private :: n_tracks ! total number of tracks + real(8), private, allocatable :: coords(:,:) ! track coordinates #endif contains !=============================================================================== ! INITIALIZE_PARTICLE_TRACK opens a particle track output file. -! -! TODO: This subroutine needs to be modified to work with HDF5 files. It -! should also probably write a header that identifies the file as a particle -! track and maybe adds particle identifying information (batch #, etc.). !=============================================================================== subroutine initialize_particle_track() character(MAX_FILE_LEN) :: filename - ! Replace the following with '#ifdef HDF5'. #if 0 - integer(HSIZE_T) :: dims(2), max_dims(2), i + integer(HSIZE_T) :: dims(2), max_dims(2), i, j(1) integer(HID_T) :: cparms - integer(HID_T) :: dapl -#endif - filename = trim(path_output) // 'track_' // trim(to_str(current_batch)) & - // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) - - ! Replace the following with '#ifdef HDF5'. -#if 0 - filename = filename // '.h5' + // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) & + // '.h5' + n_tracks = 0 + ! Create file. call h5fcreate_f(filename, H5F_ACC_TRUNC_F, track_fh, hdf5_err) + ! Create a dataspace with an unlimited max dimension. dims = (/3, 1/) i = 3 max_dims = (/i, H5S_UNLIMITED_F/) call h5screate_simple_f(2, dims, track_dspace, hdf5_err, max_dims) + ! Set dataspace chunking. call h5pcreate_f(H5P_DATASET_CREATE_F, cparms, hdf5_err) call h5pset_chunk_f(cparms, 2, dims, hdf5_err) + ! Set fill value. call h5pset_fill_value_f(cparms, H5T_NATIVE_DOUBLE, 0, hdf5_err) + ! Create dataset. call h5dcreate_f(track_fh, 'coordinates', H5T_NATIVE_DOUBLE, track_dspace, & track_dset, hdf5_err, cparms) - n_tracks = 1 +#endif +#ifdef HDF5 + filename = trim(path_output) // 'track_' // trim(to_str(current_batch)) & + // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) & + // '.h5' + n_tracks = 0 + fname = filename + allocate(coords(3,1)) #else - filename = filename // '.binary' + filename = trim(path_output) // 'track_' // trim(to_str(current_batch)) & + // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p % id)) & + // '.binary' open(UNIT=UNIT_TRACK, FILE=filename, ACTION="write", & STATUS='replace', ACCESS='stream') #endif - end subroutine initialize_particle_track !=============================================================================== ! WRITE_PARTICLE_TRACK outputs particle position to a binary file. -! -! TODO: This subroutine needs to be modified to work with HDF5 files. Perhaps -! it should also be made somehow more general so that it can output -! information other than just particle position. !=============================================================================== subroutine write_particle_track() - ! Replace the following with '#ifdef HDF5'. #if 0 integer(HSIZE_T) :: dims(2), max_dims(2), offset(2), i integer(HSIZE_T), parameter :: cnt(2) = (/3, 1/), write_dims(2) = (/3, 1/) - real(8) :: coords(3, 1) - integer(HSIZE_T) :: test + real(8) :: coords(3) + ! There is another set of track coordinates. Incriment the counter. + n_tracks = n_tracks + 1 + ! Make the dataset the right size to fit the track coordinates. i = 3 dims = (/i, n_tracks/) - call h5dget_space_f(track_dset, track_dspace, hdf5_err) call h5dset_extent_f(track_dset, dims, hdf5_err) + ! Get the dataspace with the updated dimensions. call h5dget_space_f(track_dset, track_dspace, hdf5_err) + ! Select the hyperslab where the latest coordinates will be written. i = 1 offset = (/i, n_tracks/) + if (n_tracks < 20) then call h5sselect_hyperslab_f(track_dspace, H5S_SELECT_SET_F, offset, cnt, & hdf5_err) - !call h5sget_select_npoints_f(track_dspace, test, hdf5_err) - print *, test - coords(:, 1) = p % coord0 % xyz - print *, coords - if (n_tracks<20) then + endif + ! Write the coordinates to the dataset. + coords = p % coord0 % xyz call h5dwrite_f(track_dset, H5T_NATIVE_DOUBLE, coords, & write_dims, hdf5_err, file_space_id=track_dspace) - end if - !call h5sselect_hyperslab_f(track_dspace, H5S_SELECT_SET_F, offset, dims, & - ! hdf5_err) - !i = 3 +#endif +#ifdef HDF5 + real(8) :: old_coords(3, n_tracks) + + old_coords = coords n_tracks = n_tracks + 1 - !dims = (/i, n_tracks/) - !call h5sset_extent_simple_f(track_dspace, 2, dims, max_dims, hdf5_err) - !call h5dset_extent_f(track_dset, dims, hdf5_err) + deallocate(coords) + allocate(coords(3, n_tracks)) + coords(:, 1:n_tracks-1) = old_coords + coords(:, n_tracks) = p % coord0 % xyz #else write(UNIT_TRACK) p % coord0 % xyz #endif @@ -112,24 +138,20 @@ contains !=============================================================================== ! FINALIZE_PARTICLE_TRACK closes the particle track file. -! -! TODO: This subroutine needs to be modified to work with HDF5 files. !=============================================================================== subroutine finalize_particle_track() - ! Replace the following with '#ifdef HDF5'. #if 0 - integer(HSIZE_T) :: dims(2), max_dims(2), i - - call h5sget_simple_extent_dims_f(track_dspace, dims, max_dims, hdf5_err) - !i = 3 - !dims = (/i, dims(2)-1/) - !call h5dextend_f(track_dset, dims, hdf5_err) - call h5dclose_f(track_dset, hdf5_err) call h5sclose_f(track_dspace, hdf5_err) call h5pclose_f(cparms, hdf5_err) call h5fclose_f(track_fh, hdf5_err) +#endif +#ifdef HDF5 + call file_create(fname, 'serial') + call write_data(coords, 'coordinates', length=(/3, n_tracks/)) + call file_close('serial') + deallocate(coords) #else close(UNIT=UNIT_TRACK) #endif diff --git a/src/utils/track.py b/src/utils/track.py index 1907d1ec0a..f0b1f23e5a 100755 --- a/src/utils/track.py +++ b/src/utils/track.py @@ -5,7 +5,7 @@ Usage information can be obtained by running 'track.py --help': usage: track.py [-h] [-o OUT] IN [IN ...] - Convert binary particle track file to a .pvtp file. + Convert particle track file to a .pvtp file. positional arguments: IN Input particle track data filename(s). @@ -26,9 +26,8 @@ import vtk def _parse_args(): # Create argument parser. parser = argparse.ArgumentParser( - description='Convert binary particle track file to a .pvtp file.') - parser.add_argument('input', metavar='IN', type=argparse.FileType('rb'), - nargs='+', + description='Convert particle track file to a .pvtp file.') + parser.add_argument('input', metavar='IN', type=str, nargs='+', help='Input particle track data filename(s).') parser.add_argument('-o', '-out', '--out', metavar='OUT', type=str, dest='out', @@ -39,8 +38,15 @@ def _parse_args(): def main(): - # Parse commandline argumetns. + # Parse commandline arguments. args = _parse_args() + + # Check input file extensions. + for fname in args.input: + if len(fname) > 3 and fname[-3:] == '.h5': continue + if len(fname) > 7 and fname[-7:] == '.binary': continue + raise ValueError("Input file names must either end with '.h5' or " + "'.binary'.") # Make sure that the output filename ends with '.pvtp'. if not args.out: @@ -48,19 +54,28 @@ def main(): elif os.path.splitext(args.out)[1] != '.pvtp': args.out = ''.join([args.out, '.pvtp']) + # Import HDF library if HDF files are present + for fname in args.input: + if fname[-3:] != '.h5': continue + import h5py + break + points = vtk.vtkPoints() cells = vtk.vtkCellArray() j = 0 k = 0 - for fin in args.input: - track = fin.read() - coords = [struct.unpack("ddd", track[24*i : 24*(i+1)]) - for i in range(len(track)/24)] - n_points = len(coords) - - for triplet in coords: - points.InsertNextPoint(triplet) - + for fname in args.input: + if len(fname) > 7 and fname[-7:] == '.binary': + track = open(fname, 'rb').read() + coords = [struct.unpack("ddd", track[24*i : 24*(i+1)]) + for i in range(len(track)/24)] + n_points = len(coords) + for triplet in coords: points.InsertNextPoint(triplet) + elif fname[-3:] == '.h5': + coords = h5py.File(fname).get('coordinates') + n_points = coords.shape[0] + for i in range(n_points): points.InsertNextPoint(coords[i,:]) + # Create VTK line and assign points to line. line = vtk.vtkPolyLine() line.GetPointIds().SetNumberOfIds(n_points)