Streamlines of N-body simulations with Py-SPHViewer

Figure 1) Streamlines of the velocity field of a cosmological N-body simulation. Although we can not see the density field, we can infer the location of voids and filaments from the shape of the streamlines.

Have you ever wondered how to get the streamlines of the velocity field of a N-body simulation? In practice, streamlines are an extraordinary tool to get insight on what is going on in a simulation by inspecting a single snapshot. They are also extremely useful in many applications that are outside the scope of this post.

Let me just show an example of their power. Figure 1 shows the streamlines of the velocity field of a cosmological N-body simulation. Although the underlying density field is not shown, we can infer its presence by looking at the flux lines. In fact, we can infer the presence of a “void” at (x,y) = (10,-20) because the flux lines point outward. We can also infer the location of filaments by looking at the intersection of opposing lines.

Doing this kind of diagrams is relatively straightforward once we have the velocity field calculated on a regular grid. However, for N-body simulations (or hydrodynamical SPH simulations) it is not always the case. Here is where Py-SPHViewer can give us a hand.

With Py-SPHViewer we are able to calculate the instantaneous (and continuous) mean velocity field, which can then be passed to matplotlib to construct the streamlines.

In this post, I’ll give the recipe to get incredible nice results in only a few steps.

Getting the data for this example and drawing the streamlines

In this example I use an output of a low-resolution N-Body cosmological simulation run with Gadget-2.

Figure 2) Scatter plot of the low-resolution cosmological simulation used in this example.

The initial conditions were created using Music, and I used my own laptop to run it. You may want to download this file to continue with the example. Figure 2 shows a scatter plot of the data.

In order to get the velocity field of the box we can consider different slices, or just project the full box. I’ll use the last approach, and consider the mean (projected) velocity field. In particular, I will consider the density-weighted average. Because of this, the first thing to do is to calculate the density field of the box:

This can be done with QuickView:

import h5py
from import QuickView
#Read the file
snap = h5py.File('darkmatter_box.h5py', 'r')
pos = snap['PartType1/Coordinates'].value
vel = snap['PartType1/Velocities'].value
hsml = snap['PartType1/SmoothingLength'].value
#Calculate the density field
qv = QuickView(pos.T, hsml=hsml, r='infinity', x=32, y=32, z=32,plot=False, extent=[-32,32,-32,32], logscale=False)
density_field = qv.get_image()
extent = qv.get_extent()

Fig. 3 shows the resulting density field. We now need to do the same for the x and y components of the velocity field:

Figure 3) Continuous density field of the simulated box. We will use it to calculate the density-weighted velocity field
vfield = []
for i in xrange(2):
    qv = QuickView(pos.T, vel[:,i], hsml=hsml, r='infinity', x=32, y=32, z=32,plot=False, extent=[-32,32,-32,32], logscale=False)

Where we have smoothed the velocity field (x and y components separately), and divided the resulting field by the density field, which effectively produces a density-weighted average. The resulting fields are shown in Fig. 4

Figure 4) X (left) and Y (right) components of the density-weighted average velocity field of the cosmological simulation considered in the current example. We will use this two field to calculate the corresponding streamlines.

We have now all the ingredients to let matplotlib do the magic. I will use now the streamplot method:

fig = plt.figure(1, figsize=(7,7))
ax = fig.add_subplot(111)

X = np.linspace(extent[0], extent[1], 500)
Y = np.linspace(extent[2], extent[3], 500)
ax.imshow(np.log10(density_field), origin='lower', extent=extent, cmap='bone')

v = np.log10(np.sqrt(vfield[0]**2+vfield[1]**2))
color = v/np.max(v)

lw = 2*color

streams = ax.streamplot(X,Y,vfield[0], vfield[1], color=color, density=1.5, cmap='jet', linewidth=lw, arrowsize=1)
ax.set_xlabel(r'$\rm X / \ Mpc \ h^{-1}$', size=25)
ax.set_ylabel(r'$\rm Y / \ Mpc \ h^{-1}$', size=25)

fig.savefig('figura5.png', bbox_inches='tight')

which produces the output shown in Fig. 5.


Streamlines are straight forward to calculate for N-body simulations, provided the velocity field is known on a regular grid. The velocity field can be calculated trivially using Py-SPHViewer, and streamlines can be constructed with matplotlib.

Finally, let me just say that streamlines can also be animated to give a really deep insight on what is going on in the simulation. I plan to focus on this in a forthcoming post. Please, do not hesitate to leave your comments, suggestions or questions down below.

Figure 5) Streamlines of the density-weighted average velocity field of the full box. Colours indicate the magnitude of the modulus of the velocity, from blue to red. The width of the lines are also a function of this quantity.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s