In [1]:
# Author: Volker Hoffmann <volmker@cheleb.net>, <volker@physik.uzh.ch>
# Date: 02 December 2014

In [2]:
%matplotlib inline

In [3]:
import matplotlib as mpl; mpl.rcParams['savefig.dpi'] = 144
import matplotlib.pyplot as plt
import numpy as np


In this notebook, we will use Pymses to load raw Ramses data. We generate a set of points, and ask Pymses to return the hydro fields at these locations. We discuss the alternate option of returning the hydro fields for all cells in a separate notebook. For more information, refer to the Pymses webpage and documentation at [1]. We also use the cell magic %%time below. There are whole lot more you can use [2].

In [4]:
import pymses


We need to generate the array of points we wish to get the hydro data from. We generate point coordinates as 3D array, which we then arrange into the form Pymses accepts.

In [5]:
# Generate Sampling Points. Generates 64 numbers in the range (0,1) along each dimension.
x, y, z = np.mgrid[0:1:64j, 0:1:64j, 0:1:64j]

# Reshape
npoints = np.prod(x.shape)
x1 = np.reshape(x, npoints)
y1 = np.reshape(y, npoints)
z1 = np.reshape(z, npoints)

# Arrange for Pymses
pxyz = np.array([x1, y1, z1])
pxyz = pxyz.transpose()


We now load the Ramses data, and then sample the hydro fields at the requested locations. Passing verbose=False to RamsesOutput() suppresses Pymses informing us about what files it's reading. This list can be painfully long for MPI runs.

In [6]:
%%time
# Prepare Ramses Output
basedir = "/zbox/data/volker/PymsesDemo"
output = pymses.RamsesOutput("%s" % basedir, 148, verbose=False)
source = output.amr_source(["rho", "vel", "P"])

# Sample Hydro Fields
dset = pymses.analysis.sample_points(source, pxyz, use_C_code=True)

Computing hilbert minimal domain description for output 148 ...
Done !
Warning : 5 variables found - Using default value for pymses_field_descrs.py because  [Errno 2] No such file or directory: './pymses_field_descrs.py'
CPU times: user 3.11 s, sys: 304 ms, total: 3.42 s
Wall time: 3.61 s



Pymses returns all data in a flat list. We must reshape this again for further use.

In [7]:
rho = np.reshape(dset["rho"], (64,64,64))
vx = np.reshape(dset["vel"][:,0], (64,64,64))
vy = np.reshape(dset["vel"][:,1], (64,64,64))
vz = np.reshape(dset["vel"][:,2], (64,64,64))


We compute the surface density by integrating along the z direction.

## Compute & Plot Surface Density¶

In [8]:
dz = 1.0/64.0
Sigma = np.sum(rho*dz, axis=2)

In [9]:
# Plot Surface Density
fig, ax = plt.subplots(1,1)
im = ax.imshow(np.log10(Sigma), interpolation='none', cmap="bone", extent=[-0.5,0.5,-0.5,0.5])
plt.colorbar(im)
ax.set_xlabel("X")
ax.set_ylabel("X")
ax.set_title("Surface Density")
pass


## Plotting 3D¶

3D plotting capabilities in Matplotlib are somewhat crippled. We demostrate a quiverplot for the velocity field below. If you want something more complicated (proper raytracting, contour surfaces, etc), have a look at Mayavi [1] (which at present is tricky to embed in IPython notebook).

In [10]:
from mpl_toolkits.mplot3d import Axes3D


We only consider every 8th position to now oversaturate our senses with arrows. Again, to make nice 3D plots (e.g., velocity fields superimposed on contour surfaces), have a look at Mayavi, which is built for it.

In [11]:
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.quiver(x[::8,::8,::8]-0.5, y[::8,::8,::8]-0.5, z[::8,::8,::8]-0.5, \
vx[::8,::8,::8], vy[::8,::8,::8], vz[::8,::8,::8], \
length=0.05, arrow_length_ratio=0.5, alpha=0.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Velocity Field")
pass