Merge pull request #1180 from pshriwise/zero_copy_voxel_to_vtk

Reduced Memory (and faster) Voxel to VTK

include/openmc/constants.h

@@ -33,7 +33,7 @@ constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
 constexpr std::array<int, 2> VERSION_TRACK {2, 0};
 constexpr std::array<int, 2> VERSION_SUMMARY {6, 0};
 constexpr std::array<int, 2> VERSION_VOLUME {1, 0};
-constexpr std::array<int, 2> VERSION_VOXEL {1, 0};
+constexpr std::array<int, 2> VERSION_VOXEL {2, 0};
 constexpr std::array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr char VERSION_MULTIPOLE[] {"v0.2"};
 

scripts/openmc-voxel-to-vtk

@@ -8,6 +8,8 @@ import numpy as np
 import h5py
 import vtk
 
+_min_version = (2,0)
+
 
 def main():
     # Process command line arguments
@@ -19,29 +21,40 @@ def main():
 
     # Read data from voxel file
     fh = h5py.File(args.voxel_file, 'r')
+
+    # check version
+    version = tuple(fh.attrs['version'])
+    if version < _min_version:
+        old_version = ".".join(map(str,version))
+        min_version = ".".join(map(str,_min_version))
+        err_msg = "This voxel file's version is {}. This script " \
+                  "only supports voxel files with version {} or " \
+                  "higher. Please generate a new voxel file using " \
+                  "a newer version of OpenMC.".format(old_version, min_version)
+        raise ValueError(err_msg)
+
     dimension = fh.attrs['num_voxels']
     width = fh.attrs['voxel_width']
     lower_left = fh.attrs['lower_left']
-    voxel_data = fh['data'].value
 
     nx, ny, nz = dimension
     upper_right = lower_left + width*dimension
-    
+
     grid = vtk.vtkImageData()
     grid.SetDimensions(nx+1, ny+1, nz+1)
     grid.SetOrigin(*lower_left)
     grid.SetSpacing(*width)
-    
-    data = vtk.vtkDoubleArray()
+
+    # transpose data from OpenMC ordering (zyx) to VTK ordering (xyz)
+    # and flatten to 1-D array
+    print("Reading and translating data...")
+    h5data = fh['data'][...]
+
+    data = vtk.vtkIntArray()
     data.SetName("id")
-    data.SetNumberOfTuples(nx*ny*nz)
-    for x in range(nx):
-        sys.stdout.write(" {}%\r".format(int(x/nx*100)))
-        sys.stdout.flush()
-        for y in range(ny):
-            for z in range(nz):
-                i = z*nx*ny + y*nx + x
-                data.SetValue(i, voxel_data[x, y, z])
+    # set the array using the h5data array
+    data.SetArray(h5data, h5data.size, True)
+    # add data to image grid
     grid.GetCellData().AddArray(data)
 
     writer = vtk.vtkXMLImageDataWriter()
@@ -52,6 +65,7 @@ def main():
     if not args.output.endswith(".vti"):
         args.output += ".vti"
     writer.SetFileName(args.output)
+    print("Writing VTK file {}...".format(args.output))
     writer.Write()
 
 if __name__ == '__main__':

src/plot.cpp

@@ -884,9 +884,9 @@ void create_voxel(Plot pl)
   // Create dataset for voxel data -- note that the dimensions are reversed
   // since we want the order in the file to be z, y, x
   hsize_t dims[3];
-  dims[0] = pl.pixels_[0];
+  dims[0] = pl.pixels_[2];
   dims[1] = pl.pixels_[1];
-  dims[2] = pl.pixels_[2];
+  dims[2] = pl.pixels_[0];
   hid_t dspace, dset, memspace;
   voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
 
@@ -895,33 +895,33 @@ void create_voxel(Plot pl)
   ll[1] = ll[1] + vox[1] / 2.;
   ll[2] = ll[2] + vox[2] / 2.;
 
-  int data[pl.pixels_[1]][pl.pixels_[2]];
+  int data[pl.pixels_[1]][pl.pixels_[0]];
 
   ProgressBar pb;
 
   RGBColor rgb;
   int id;
-  for (int x = 0; x < pl.pixels_[0]; x++) {
-    pb.set_value(100.*(double)x/(double)(pl.pixels_[0]-1));
+  for (int z = 0; z < pl.pixels_[2]; z++) {
+    pb.set_value(100.*(double)z/(double)(pl.pixels_[2]-1));
     for (int y = 0; y < pl.pixels_[1]; y++) {
-      for (int z = 0; z < pl.pixels_[2]; z++) {
+      for (int x = 0; x < pl.pixels_[0]; x++) {
         // get voxel color
         position_rgb(p, pl, rgb, id);
         // write to plot data
-        data[y][z] = id;
-        // advance particle in z direction
-        p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2];
+        data[y][x] = id;
+        // advance particle in x direction
+        p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0];
       }
       // advance particle in y direction
       p.coord[0].xyz[1] = p.coord[0].xyz[1] + vox[1];
-      p.coord[0].xyz[2] = ll[2];
+      p.coord[0].xyz[0] = ll[0];
     }
-    // advance particle in x direction
-    p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0];
+    // advance particle in z direction
+    p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2];
     p.coord[0].xyz[1] = ll[1];
-    p.coord[0].xyz[2] = ll[2];
+    p.coord[0].xyz[0] = ll[0];
     // Write to HDF5 dataset
-    voxel_write_slice(x, dspace, dset, memspace, &(data[0]));
+    voxel_write_slice(z, dspace, dset, memspace, &(data[0]));
   }
 
   voxel_finalize(dspace, dset, memspace);

tests/regression_tests/plot/results_true.dat

@@ -1 +1 @@
-01ecda0f3820a49c8a41d8dc47d1e5c58767a04301621c2437231fcc04401ddea47b67d0529ca56a32d4d97b4f1416a2e0b6120d3bdc87d74a7e9889758a8808
+30f435795c755b6fb28565a33c68e95415e5c09bca59f5549a6cbd2dd165bac9d90f5a599e12a5041cb927698abde562b7bac14175441c57b5c0dccd00544506

tests/regression_tests/plot_voxel/results_true.dat

@@ -1 +1 @@
-20a0c3598bb698efa4bba65c5e55d71f4385e5b1bcf0067964e45001c12f5c0a729935a96409c760dbd3ec439ae6c676fce77d6e578b41f8f3e0ac68c9277acb
+76274390783bf0fd54a63e3ee141e072e797a564935a16f22a97788bc051e33fd7160dd311f0c50e45b0b498701f6735a37fabacbf9f36903381d9e18759a684