mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-25 20:45:35 -04:00
Use HDF5 format for voxel file rather than raw binary
This commit is contained in:
Paul Romano
parent
8b67fa7a92
commit
3274cf3f6d
4 changed files with 72 additions and 55 deletions
|
|
@ -12,3 +12,4 @@ Output File Formats
|
|||
source
|
||||
particle_restart
|
||||
track
|
||||
voxel
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
54
src/plot.F90
54
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
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue