Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Replacing legacy vtk file writing with h5py file writing for meshes, geometry and tracks for paraview > 5.13.0 #3113

Open
shimwell opened this issue Aug 9, 2024 · 2 comments

Comments

@shimwell
Copy link
Member

shimwell commented Aug 9, 2024

Description

I'm attending at a Paraview presentation by @mwestphal and heard about this nice new feature that I think we should make use of in openmc.

Paraview is supporting hdf5 files if written in a specific way with h5py 🎉 .

Here is the relevant blog post https://www.kitware.com/how-to-write-time-dependent-data-in-vtkhdf-files/ with a nice example of writing the data using h5py

This is interesting for me because:

  • We currently have a vtk dependency which could be removed
  • We already have h5py as a dependency which can be utilized more
  • The vtk files we currently write are legacy format and not recommend by paraview anymore
  • The hdf5 files we could write for our 3d data to would be smaller in size and quicker to load than the .vtk files

The parts of the code that could be updated are

Alternatives

Stick with the current solutions, or add an arg to allow "legacy" vtk file writing

Compatibility

Will require ParaView 5.13.0 or newer to open the hd5f based files

@shimwell
Copy link
Member Author

shimwell commented Jan 7, 2025

I found a few examples for writing VTKHDF files here

I adapted one of the examples so that it writes random data to a mesh

Still got to adapt this so that it takes points, coordinates etc from the mesh but I think this is a reasonable start

import numpy as np
import h5py as h5


# Convenient method to fill h5py node with a numpy array
def append_dataset(dset, array):
    origLen = dset.shape[0]
    dset.resize(origLen + array.shape[0], axis=0)
    dset[origLen:] = array
    if len(array) == 1:
        print(f"Appended {len(array)} items {(array)} to {dset.name}, new length: {dset.shape[0]}")
    else:
        print(f"Appended {len(array)} items to {dset.name}, new length: {dset.shape[0]}")


# Fill VTKHDF file with cubes
def fill_with_dummy_unstructured_grid(root, nCubePerDim = 1):
    # Additional meta information to generate a valid VTKHDF file
    fullDim = nCubePerDim * nCubePerDim * nCubePerDim
    cubePoints = 8
    numberOfConnectivity = 8
    numberOfPts = fullDim * cubePoints
    numberOfOffset = fullDim + 1

    points = np.array(np.empty([numberOfPts,3]))
    connectivity = np.array(np.empty([numberOfConnectivity * fullDim]))
    offsets = np.array(np.empty([fullDim+1]))
    types = np.array(np.empty([numberOfOffset-1]))

    [fillGeometry(idx, points, connectivity, offsets, types, numberOfConnectivity, nCubePerDim) for idx in range(0,fullDim)]

    offsets[fullDim] = fullDim * 8
    # print('offsets', offsets)


    append_dataset(root['NumberOfPoints'], np.array([numberOfPts]))
    append_dataset(root['Points'], points)

    append_dataset(root['NumberOfConnectivityIds'], np.array([numberOfConnectivity*fullDim]))
    append_dataset(root['Connectivity'], connectivity)

    append_dataset(root['NumberOfCells'], np.array([numberOfOffset-1]))
    append_dataset(root['Offsets'], offsets)

    append_dataset(root['Types'], types)

    # Create and populate the scalar data values dataset
    scalar_data = np.random.rand(numberOfOffset-1)  # Example: random scalar values
    cell_data_group = root.create_group('CellData')

    cell_data_group.create_dataset('Scalars', (0,), maxshape=(None,), dtype='float64', chunks=True)

    append_dataset(cell_data_group['Scalars'], scalar_data)


def fillGeometry(cubeIdx, points, connectivity,offsets, types, numberOfConnectivity, nCubePerDim):

    XIdx = cubeIdx % nCubePerDim
    XQuot = cubeIdx // nCubePerDim

    YIdx = XQuot % nCubePerDim
    YQuot = XQuot // nCubePerDim

    ZIdx = YQuot % nCubePerDim

    points[0 + cubeIdx * 8] = [2 + (XIdx + 0) * 0.25 ,(YIdx + 0) * 0.25, (ZIdx + 0) * 0.25]
    points[1 + cubeIdx * 8] = [2 + (XIdx + 1) * 0.25 ,(YIdx + 0) * 0.25, (ZIdx + 0) * 0.25]
    points[2 + cubeIdx * 8] = [2 + (XIdx + 0) * 0.25 ,(YIdx + 1) * 0.25, (ZIdx + 0) * 0.25]
    points[3 + cubeIdx * 8] = [2 + (XIdx + 1) * 0.25 ,(YIdx + 1) * 0.25, (ZIdx + 0) * 0.25]
    points[4 + cubeIdx * 8] = [2 + (XIdx + 0) * 0.25 ,(YIdx + 0) * 0.25, (ZIdx + 1) * 0.25]
    points[5 + cubeIdx * 8] = [2 + (XIdx + 1) * 0.25 ,(YIdx + 0) * 0.25, (ZIdx + 1) * 0.25]
    points[6 + cubeIdx * 8] = [2 + (XIdx + 0) * 0.25 ,(YIdx + 1) * 0.25, (ZIdx + 1) * 0.25]
    points[7 + cubeIdx * 8] = [2 + (XIdx + 1) * 0.25 ,(YIdx + 1) * 0.25, (ZIdx + 1) * 0.25]


    connectivity[0 + cubeIdx * 8] = 0 + numberOfConnectivity * cubeIdx
    connectivity[1 + cubeIdx * 8] = 1 + numberOfConnectivity * cubeIdx
    connectivity[2 + cubeIdx * 8] = 3 + numberOfConnectivity * cubeIdx
    connectivity[3 + cubeIdx * 8] = 2 + numberOfConnectivity * cubeIdx
    connectivity[4 + cubeIdx * 8] = 4 + numberOfConnectivity * cubeIdx
    connectivity[5 + cubeIdx * 8] = 5 + numberOfConnectivity * cubeIdx
    connectivity[6 + cubeIdx * 8] = 7 + numberOfConnectivity * cubeIdx
    connectivity[7 + cubeIdx * 8] = 6 + numberOfConnectivity * cubeIdx

    cubePoints = 8
    offsets[cubeIdx] = cubeIdx * cubePoints

    # Code recognized by VTK as HEXAHEDRON
    types[cubeIdx] = 12



f = h5.File('unstructured.hdf', 'w')

root = f.create_group('VTKHDF')
root.attrs['Version'] = (2,1)
ascii_type = 'UnstructuredGrid'.encode('ascii')
root.attrs.create('Type', ascii_type, dtype=h5.string_dtype('ascii', len(ascii_type)))

root.create_dataset('NumberOfPoints', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('Types', (0,), maxshape=(None,), dtype='uint8')
root.create_dataset('Points', (0,3), maxshape=(None,3), dtype='f')

root.create_dataset('NumberOfConnectivityIds', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('NumberOfCells', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('Offsets', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('Connectivity', (0,), maxshape=(None,), dtype='i8')

fill_with_dummy_unstructured_grid(root, 4)

@shimwell
Copy link
Member Author

shimwell commented Jan 7, 2025

Ok I've done a bit more work on this and I can now write out a DAGMC volume mesh using the VTKHDF format and hyp5 (no use of vtk)

I've got a branch with a working version over here develop...shimwell:openmc:adding_write_data_to_vtk_hdf

Screenshot From 2025-01-07 17-26-25

import openmc
import numpy as np
import h5py as h5

def append_dataset(dset, array):
    origLen = dset.shape[0]
    dset.resize(origLen + array.shape[0], axis=0)
    dset[origLen:] = array
    if len(array) == 1:
        print(f"Appended {len(array)} items {(array)} to {dset.name}, new length: {dset.shape[0]}")
    else:
        print(f"Appended {len(array)} items to {dset.name}, new length: {dset.shape[0]}")

def remove_trailing_negatives(connectivity):
    """
    Remove trailing -1 values from each sub-array in the connectivity array.

    Parameters:
    connectivity (np.ndarray): The connectivity array with potential trailing -1 values.

    Returns:
    np.ndarray: The cleaned connectivity array with trailing -1 values removed.
    """
    cleaned_connectivity = []
    for cell in connectivity:
        # Find the index of the first -1 value, if any
        first_negative_index = np.where(cell == -1)[0]
        if first_negative_index.size > 0:
            # Slice the array up to the first -1 value
            cleaned_connectivity.append(cell[:first_negative_index[0]])
        else:
            # No -1 values, append the whole cell
            cleaned_connectivity.append(cell)
    
    return np.array(cleaned_connectivity, dtype='int32')

sp=openmc.StatePoint('statepoint.100.h5')

um = sp.meshes[1]  # might need to change this index depending on the mesh you want to convert

trimmed_connectivity=remove_trailing_negatives(um.connectivity)
connectivity = trimmed_connectivity.flatten()

points_per_cell=4

offsets = np.arange(0, um.n_elements * points_per_cell+1, points_per_cell)

f = h5.File('unstructured_from_openmc.hdf', 'w')


root = f.create_group('VTKHDF')
root.attrs['Version'] = (2,1)
ascii_type = 'UnstructuredGrid'.encode('ascii')
root.attrs.create('Type', ascii_type, dtype=h5.string_dtype('ascii', len(ascii_type)))

root.create_dataset('NumberOfPoints', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('Types', (0,), maxshape=(None,), dtype='uint8')
root.create_dataset('Points', (0,3), maxshape=(None,3), dtype='f')
root.create_dataset('NumberOfConnectivityIds', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('NumberOfCells', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('Offsets', (0,), maxshape=(None,), dtype='i8')
root.create_dataset('Connectivity', (0,), maxshape=(None,), dtype='i8')

append_dataset(root['NumberOfPoints'], np.array([len(um.vertices)]))
append_dataset(root['Points'], um.vertices)

append_dataset(root['NumberOfConnectivityIds'], np.array([len(connectivity)]))
append_dataset(root['Connectivity'], connectivity)

append_dataset(root['NumberOfCells'], np.array([um.n_elements]))
append_dataset(root['Offsets'], offsets)

types = np.full(um.n_elements, 10, dtype='uint8')  # VTK_TETRA type assumes DAGMC mesh
append_dataset(root['Types'], types)

scalar_data = np.random.rand(um.n_elements)  # Example: random scalar values
cell_data_group = root.create_group('CellData')

cell_data_group.create_dataset('Scalars', (0,), maxshape=(None,), dtype='float64', chunks=True)

append_dataset(cell_data_group['Scalars'], scalar_data)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant