From 3274cf3f6dcdccfa5dd55135f699bfeb33cbe03a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 17 Sep 2015 12:47:20 +0700 Subject: [PATCH] Use HDF5 format for voxel file rather than raw binary --- docs/source/usersguide/output/index.rst | 1 + scripts/openmc-voxel-to-silovtk | 69 +++++++++---------------- src/constants.F90 | 3 +- src/plot.F90 | 54 +++++++++++++++---- 4 files changed, 72 insertions(+), 55 deletions(-) diff --git a/docs/source/usersguide/output/index.rst b/docs/source/usersguide/output/index.rst index 1eb85e9d5..161a8de2d 100644 --- a/docs/source/usersguide/output/index.rst +++ b/docs/source/usersguide/output/index.rst @@ -12,3 +12,4 @@ Output File Formats source particle_restart track + voxel diff --git a/scripts/openmc-voxel-to-silovtk b/scripts/openmc-voxel-to-silovtk index eb052c75f..1329b6b6c 100755 --- a/scripts/openmc-voxel-to-silovtk +++ b/scripts/openmc-voxel-to-silovtk @@ -1,9 +1,11 @@ -#!/usr/bin/env python2 +#!/usr/bin/env python from __future__ import division, print_function import struct import sys +import numpy as np +import h5py def parse_options(): """Process command line arguments""" @@ -22,14 +24,16 @@ def parse_options(): return parsed -def main(file_, o): - print(file_) - fh = open(file_, 'rb') - header = get_header(fh) - meshparms = (header['dimension'] + header['lower_left'] + - header['upper_right']) - nx, ny, nz = meshparms[:3] - ll = header['lower_left'] +def main(filename, o): + # Read data from voxel file + fh = h5py.File(filename, 'r') + dimension = fh['num_voxels'].value + width = fh['voxel_width'].value + lower_left = fh['lower_left'].value + voxel_data = fh['data'].value + + nx, ny, nz = dimension + upper_right = lower_left + width*dimension if o.vtk: try: @@ -40,13 +44,10 @@ def main(file_, o): 'See: http://www.vtk.org/') return - origin = [(l + w*n/2.) for n, l, w in - zip((nx, ny, nz), ll, header['width'])] - grid = vtk.vtkImageData() grid.SetDimensions(nx+1, ny+1, nz+1) - grid.SetOrigin(*ll) - grid.SetSpacing(*header['width']) + grid.SetOrigin(*lower_left) + grid.SetSpacing(*width) data = vtk.vtkDoubleArray() data.SetName("id") @@ -57,8 +58,7 @@ def main(file_, o): for y in range(ny): for z in range(nz): i = z*nx*ny + y*nx + x - id_ = get_int(fh)[0] - data.SetValue(i, id_) + data.SetValue(i, voxel_data[x,y,z]) grid.GetCellData().AddArray(data) writer = vtk.vtkXMLImageDataWriter() @@ -81,44 +81,23 @@ def main(file_, o): if not o.output.endswith(".silo"): o.output += ".silo" silomesh.init_silo(o.output) - silomesh.init_mesh('plot', *meshparms) + meshparams = list(map(int, dimension)) + list(map(float, lower_left)) + \ + list(map(float, upper_right)) + silomesh.init_mesh('plot', *meshparams) silomesh.init_var("id") - for x in range(1, nx+1): + for x in range(nx): sys.stdout.write(" {0}%\r".format(int(x/nx*100))) sys.stdout.flush() - for y in range(1, ny+1): - for z in range(1, nz+1): - id_ = get_int(fh)[0] - silomesh.set_value(float(id_), x, y, z) + for y in range(ny): + for z in range(nz): + silomesh.set_value(float(voxel_data[x,y,z]), + x + 1, y + 1, z + 1) print() silomesh.finalize_var() silomesh.finalize_mesh() silomesh.finalize_silo() -def get_header(file_): - nx, ny, nz = get_int(file_, 3) - wx, wy, wz = get_double(file_, 3) - lx, ly, lz = get_double(file_, 3) - header = {'dimension': [nx, ny, nz], 'width': [wx, wy, wz], - 'lower_left': [lx, ly, lz], - 'upper_right': [lx+wx*nx, ly+wy*ny, lz+wz*nz]} - return header - - -def get_data(file_, n, typeCode, size): - return list(struct.unpack('={0}{1}'.format(n, typeCode), - file_.read(n*size))) - - -def get_int(file_, n=1, path=None): - return get_data(file_, n, 'i', 4) - - -def get_double(file_, n=1, path=None): - return get_data(file_, n, 'd', 8) - - if __name__ == '__main__': (options, args) = parse_options() if args: diff --git a/src/constants.F90 b/src/constants.F90 index 6be1326a8..3d7ca890d 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -20,7 +20,8 @@ module constants FILETYPE_STATEPOINT = -1, & FILETYPE_PARTICLE_RESTART = -2, & FILETYPE_SOURCE = -3, & - FILETYPE_TRACK = -4 + FILETYPE_TRACK = -4, & + FILETYPE_VOXEL = -5 ! ============================================================================ ! ADJUSTABLE PARAMETERS diff --git a/src/plot.F90 b/src/plot.F90 index cd4e8642e..e5b265c5c 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -5,6 +5,7 @@ module plot use geometry, only: find_cell, check_cell_overlap use geometry_header, only: Cell, BASE_UNIVERSE use global + use hdf5_interface use mesh, only: get_mesh_indices use output, only: write_message use particle_header, only: Particle, LocalCoord @@ -14,6 +15,8 @@ module plot use progress_header, only: ProgressBar use string, only: to_str + use hdf5 + implicit none contains @@ -347,11 +350,20 @@ contains integer :: x, y, z ! voxel location indices integer :: rgb(3) ! colors (red, green, blue) from 0-255 integer :: id ! id of cell or material - integer :: unit_plot ! voxel file unit + integer :: hdf5_err + integer, target :: data(pl%pixels(3),pl%pixels(2)) + integer(HID_T) :: file_id + integer(HID_T) :: dspace + integeR(HID_T) :: memspace + integer(HID_T) :: dset + integer(HSIZE_T) :: dims(3) + integer(HSIZE_T) :: dims_slab(3) + integer(HSIZE_T) :: offset(3) real(8) :: vox(3) ! x, y, and z voxel widths real(8) :: ll(3) ! lower left starting point for each sweep direction type(Particle) :: p type(ProgressBar) :: progress + type(c_ptr) :: f_ptr ! compute voxel widths in each direction vox = pl % width/dble(pl % pixels) @@ -366,11 +378,30 @@ contains p % coord(1) % universe = BASE_UNIVERSE ! Open binary plot file for writing - open(NEWUNIT=unit_plot, FILE=pl % path_plot, STATUS='replace', & - ACCESS='stream') + file_id = file_create(pl%path_plot) ! write plot header info - write(unit_plot) pl % pixels, vox, ll + call write_dataset(file_id, "filetype", FILETYPE_VOXEL) + call write_dataset(file_id, "num_voxels", pl%pixels) + call write_dataset(file_id, "voxel_width", vox) + call write_dataset(file_id, "lower_left", ll) + + ! Create dataset for voxel data -- note that the dimensions are reversed + ! since we want the order in the file to be z, y, x + dims(:) = [pl%pixels(3), pl%pixels(2), pl%pixels(1)] + call h5screate_simple_f(3, dims, dspace, hdf5_err) + call h5dcreate_f(file_id, "data", H5T_NATIVE_INTEGER, dspace, dset, hdf5_err) + + ! Create another dataspace for 2D array in memory + dims_slab(1) = pl%pixels(3) + dims_slab(2) = pl%pixels(2) + dims_slab(3) = 1 + call h5screate_simple_f(2, dims_slab(1:2), memspace, hdf5_err) + + ! Initialize offset and get pointer to data + offset(:) = 0 + call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims_slab, hdf5_err) + f_ptr = c_loc(data) ! move to center of voxels ll = ll + vox / TWO @@ -379,22 +410,19 @@ contains call progress % set_value(dble(x)/dble(pl % pixels(1))*100) do y = 1, pl % pixels(2) do z = 1, pl % pixels(3) - ! get voxel color call position_rgb(p, pl, rgb, id) ! write to plot file - write(unit_plot) id + data(z,y) = id ! advance particle in z direction p % coord(1) % xyz(3) = p % coord(1) % xyz(3) + vox(3) - end do ! advance particle in y direction p % coord(1) % xyz(2) = p % coord(1) % xyz(2) + vox(2) p % coord(1) % xyz(3) = ll(3) - end do ! advance particle in y direction @@ -402,9 +430,17 @@ contains p % coord(1) % xyz(2) = ll(2) p % coord(1) % xyz(3) = ll(3) + ! Write to HDF5 dataset + offset(3) = x - 1 + call h5soffset_simple_f(dspace, offset, hdf5_err) + call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, & + mem_space_id=memspace, file_space_id=dspace) end do - close(unit_plot) + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + call h5sclose_f(memspace, hdf5_err) + call file_close(file_id) end subroutine create_3d_dump