Basic usage of Py-SPHViewer 2.0

Introduction

This tutorial walks you through a Python-based workflow for visualizing cosmological simulation data, specifically using the Py-SPHViewer2 package for rendering smoothed particle hydrodynamics (SPH) data. We’ll also utilize h5py for handling HDF5 files.

The code reads simulation data from an HDF5 file, processes it, and then visualizes a 2D projection of the gas density field. This kind of visualization is common when looking at cosmological simulations.

Prerequisites

Before we dive into the tutorial, make sure you have the following libraries installed:

• numpy
• matplotlib
• h5py
• sphviewer2 (part of the smeshl package)

You can install these using pip:

pip install numpy matplotlib h5py py-sphviewer2

Ensure you have installed the Py-SPHViewer2 wrapper and the Smooth Mesh Library.

Step-by-Step Code Walkthrough

1. Inline Plotting and Imports

import numpy as np
import matplotlib.pyplot as plt
import h5py
from sphviewer2 import sphviewer2

• Imports:
  • h5py: This library is used to read and write HDF5 files, a common format for large numerical datasets like cosmological simulations.
  • sphviewer2: This is the Smooth Mesh Library wrapper (Py-SPHViewer 2.0) responsible for rendering SPH data into a 2D image.

2. Initialize SPH Viewer

sph = sphviewer2()

Here, we instantiate sphviewer2. If you encounter an error at this step, it is likely due to a misconfiguration of the path to the Smoothed Shared Library. Ensure that the library is either located in a standard directory or that the LD_LIBRARY_PATH (Ubuntu) or DYLD_LIBRARY_PATH (MacOS) environment variable is correctly defined.

3. Loading the Simulation Data

For this tutorial, we will use a simulation snapshot from the CAMELS project. Before continuing, download the following file:

https://users.flatironinstitute.org/~camels/Sims/IllustrisTNG/CV/CV_0/snapshot_090.hdf5

snap = h5py.File("./snapshot_090.hdf5", 'r')

Here, we load an HDF5 file that contains the snapshot of a cosmological simulation. The data might include positions, velocities, masses, and other properties of particles like gas, stars, and dark matter.

• snap: This is a file object representing the loaded HDF5 file.
• 'r': Opens the file in read-only mode, since we are only reading data.

Note: Change the file path to where your snapshot file is located.

4. Extracting Particle Data

Next, we extract specific data from the HDF5 file related to gas particles:

pgas = snap['PartType0/Coordinates'][()]
rho  = snap['PartType0/Density'][()]*1e10
mgas = snap['PartType0/Masses'][()]*1e10

• pgas: The coordinates (position) of the gas particles (PartType0 refers to gas particles in these simulations, while other types refer to dark matter, stars, etc.).
• rho: The density of the gas particles, scaled by 1e10 for unit conversion.
• mgas: The mass of the gas particles, also scaled by 1e10.

5. Calculating the Smoothing Length

hgas = (32 * mgas / (4./3.*np.pi*rho))**(1./3.)

• hgas: This is the smoothing length of the gas particles. In SPH, the smoothing length is a parameter that defines how far the particle’s influence extends. It can be estimated using the particle mass (mgas) and density (rho), assuming that the particle occupies a volume proportional to its mass and density.

The formula used is:

\(h = \left( \frac{N_{\rm ngb} m}{\frac{4}{3} \pi \rho} \right)^{\frac{1}{3}}\),

where \(N_{\rm ngb}\) is the typical number of neighbours, \(m\) is the gas particle mass, and \(\rho\) is the gas density at the particle location. In the previous example we assumed \(N_{\rm ngb}=32\).

6. Extracting the Simulation Box Size

BoxSize = snap['Header'].attrs['BoxSize']

We retrieve the size of the simulation box from the snapshot’s header. The BoxSize defines the physical dimensions of the simulation’s cubic periodic volume.

7. Rendering the Image with Py-SPHViewer 2.0

xmin = 0.0
ymin = 0.0
LevelMin = 8
LevelMax = 11
image = sph.render(pgas[:,0], pgas[:,1], mgas, hgas, xmin, ymin, BoxSize, LevelMin, LevelMax, k=None)

Here, we call the render method of the sphviewer2 object, passing in the particle data and other parameters to generate a 2D projection of the density field.

Parameters:

• pgas[:,0], pgas[:,1]: The positions of the gas particles in the x and y dimensions.
• mgas: The masses of the particles.
• hgas: The smoothing length of the particles.
• xmin, ymax: These are the minimum x and y coordinates for the image’s bounding box (bottom-left corner).
• BoxSize: The size of the simulation box, setting the upper limit of the image’s bounding box (top-right corner).
• LevelMin, LevelMax: These numbers represent the minimum and maximum levels of refinement for the adaptive mesh used in the rendering process. The higher the maximum refinement level, the better the resolution of the resulting image. The lower the minimum level, the faster the calculation.
• k=None: If needed, we can pass in the indices of a subset of particles.

The function returns an image that contains the smoothed projected density field traced by the selected SPH particles.

10. Displaying the Image

Finally, we use imshow from Matplotlib to display the rendered image. Since the density values can span many orders of magnitude, we apply a logarithmic scale (np.log10) to improve the contrast between low and high-density regions.

plt.imshow(np.log10(image))
plt.show()

The full code for this tutorial is thus:

import numpy as np
import matplotlib.pyplot as plt
import h5py
from sphviewer2 import sphviewer2

# Initialize SPH Viewer
sph = sphviewer2()

# Load the simulation data from an HDF5 file
snap = h5py.File("./snapshot_090.hdf5", 'r')

# Extract gas particle data (coordinates, density, and mass)
pgas = snap['PartType0/Coordinates'][()]
rho  = snap['PartType0/Density'][()] * 1e10
mgas = snap['PartType0/Masses'][()] * 1e10

# Calculate the smoothing length
hgas = (32 * mgas / (4./3. * np.pi * rho)) ** (1./3.)

# Extract the box size from the header of the snapshot
BoxSize = snap['Header'].attrs['BoxSize']

# Rendering parameters
xmin = 0.0
ymin = 0.0
LevelMin = 8
LevelMax = 11

# Render a 2D density projection image using Py-SPHViewer 2.0
image = sph.render(pgas[:,0], pgas[:,1], mgas, hgas, xmin, ymin, BoxSize, LevelMin, LevelMax, k=None)

# Display the rendered image with a logarithmic scale for better contrast
plt.imshow(np.log10(image))
plt.show()

and produces the following image:

Note

The beta version of the package restricts usage to levels between 8 and 11, effectively capping the image resolution at 2048 pixels or approximately 4.2 megapixels. While this resolution is adequate for most applications, obtaining the true density field traced by particles for scientific purposes requires significantly higher resolution. If this is necessary for your work, please contact me to discuss potential collaboration.

Conclusion

In this tutorial, you learned how to:

1. Load cosmological simulation data from an HDF5 file using h5py.
2. Extract particle properties like position, mass, and density.
3. Calculate the SPH smoothing length for the particles.
4. Randomly sample particles for visualization.
5. Render a 2D projection of the gas density field using Py-SPHViewer 2.0.
6. Visualize the density map with Matplotlib.

This process can be adapted to different kinds of particle data, such as dark matter or stars, and can be extended to visualize other aspects of cosmological simulations.

If you are interested in finding out how to use the ipyparallel library to speed up your calculations, check the following link:

Parallel use of Py-SPHViewer 2.0