# Author: Volker Hoffmann <firstname.lastname@example.org>, <email@example.com> # Date: 02 December 2014
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 . We also use the cell magic %%time below. There are whole lot more you can use .
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.
# 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.
%%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.
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.
dz = 1.0/64.0 Sigma = np.sum(rho*dz, axis=2)
# 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
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  (which at present is tricky to embed in IPython notebook).
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.
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