Core IO H5py

Visualizing and postprocessing particle simulation output

Contributed by: @dineshadepu

Once particle data is output from Cabana with the HDF5 interface, it can be postprocessed or visualized with many tools. This section demonstrates how to input the results, process the data, and plot an example (in this case, the position of a particle) over time using Python.

As an example, we consider the HDF5 output example, which can be run from wherever the examples were built:

./build/example/core_tutorial/13_hdf5_output/HDF5Output

This will produce a series of HDF5 (particle data) and XMF (metadata) files representing the simulation state at each output step.

Python setup

To follow this guide, you need the following:

• Python 3.10 or higher
• h5py for handling HDF5 files
• matplotlib for plotting

If needed, you can install the required libraries using pip:

pip install h5py matplotlib

Postprocessing

This analysis can be done from a script or interactively. First, load the required packages:

import h5py
import matplotlib.pyplot as plt
import os

Next, locate the output files and sort them in ascending order based on the timestep:

directory_name = "./"
files = [filename for filename in os.listdir(directory_name)
         if filename.startswith("particles") and filename.endswith("h5")]
files.sort()

files_num = []
for f in files:
    f_last = f[10:]  # Extract the step number from the filename
    files_num.append(int(f_last[:-3]))  # Remove the ".h5" extension
files_num.sort()

sorted_files = []
for num in files_num:
    sorted_files.append(f"particles_{num}.h5")

files = sorted_files

We loop through each file, extract the position data of the first particle, and plot it against time. In a more realistic case, an average or other reduction over particles would likely be of interest.

position = []
time = []
for f in files:
    with h5py.File(directory_name + f, "r") as file:
        # Access the x-position of the first particle.
        # This uses the name from the Cabana particle slice
        position.append(file["positions"][0][0])
        # Extract the simulation time from attributes
        time.append(file.attrs["Time"])

plt.plot(time, position, label="Position")
plt.xlabel("Time")
plt.ylabel("Position")
plt.legend()
plt.grid(True)
plt.show()