In [1]:
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib import cm
import numpy as np
import moments
from ppmpy import ppm
import logging
mpl_logger = logging.getLogger('matplotlib')
mpl_logger.setLevel(logging.WARNING)
moms = {}
yp = {}
for rid in ['I5', 'I11']:
moms[rid] = moments.Moments('/data/ASDR/PPMstar/C-ingestion/{:s}/FVandMoms48/'.\
format(rid), use_e3d=False)
yp[rid] = ppm.yprofile('/data/ASDR/PPMstar/C-ingestion/{:s}/YProfile/'.\
format(rid), filename_offset=-1)
In [6]:
def stretch_colormap(cmap, gamma=1., midpoint=0., N=256):
t = np.linspace(0., 1., N)
t[t < midpoint] = midpoint*(t[t < midpoint]/midpoint)**(1./gamma)
t[t > midpoint] = midpoint + (1. - midpoint)*\
((t[t > midpoint] - midpoint)/(1. - midpoint))**gamma
clrs = [cmap(tt) for tt in t]
stretched_cmap = colors.ListedColormap(clrs, 256)
return stretched_cmap
def get_Mach(moms, dump):
gamma_gas = 5./3.
p = moms.get('Prs', dump)
rho = moms.get('Rho', dump)
rhoux = moms.get('RhoUx', dump)
rhouy = moms.get('RhoUy', dump)
rhouz = moms.get('RhoUz', dump)
ux = rhoux/rho
uy = rhouy/rho
uz = rhouz/rho
nx, ny, nz = p.shape
x = np.arange(nx) - (nx - 1.)/2.
x = np.reshape(x, (nx, 1, 1))
x = np.tile(x, (1, ny, nz))
y = np.arange(ny) - (ny - 1.)/2.
y = np.reshape(y, (1, ny, 1))
y = np.tile(y, (nx, 1, nz))
z = np.arange(nz) - (nz - 1.)/2.
z = np.reshape(z, (1, 1, nz))
z = np.tile(z, (nx, ny, 1))
r = np.sqrt(x**2 + y**2 + z**2)
x /= r
y /= r
z /= r
ur = ux*x + uy*y + uz*z
utx = ux - ur*x
uty = uy - ur*y
utz = uz - ur*z
ut = np.sqrt(utx**2 + uty**2 + utz**2)
cs = np.sqrt(gamma_gas*p/rho)
Mar = np.abs(ur)/cs
Mat = np.abs(ut)/cs
return Mar, Mat
def get_rho_fluct(moms, dump):
rho = moms.get('Rho', dump)
rho_1D = moments.radprof(rho)
rho_avg = moms.fromradprof(moms.raxis, rho_1D)
drho = rho - rho_avg
# This should be a vector of zeroes, but it is not because small
# differences in the definition of the radial axis in fromradprof()
# matter. We have to do an extra iteration.
drho_1D = moments.radprof(drho)
drho_avg = moms.fromradprof(moms.raxis, drho_1D)
drho -= drho_avg
rho_fluct = drho/rho_avg
return np.array(rho_fluct)
def plot_panels(moms, dump, slice_width, vars, vlims, cmaps, axis=0, interpolation='spline16', \
ifig=1, dpi=100, cbars=True, cbar_labels=None, labels=None, label_coords=None, \
label_colors=None):
data = {}
data['Mar'], data['Mat'] = get_Mach(moms, dump)
data['rho_fluct'] = get_rho_fluct(moms, dump)
nvars = len(vars)
panel_width = 2.5
panel_height = 3.1 if cbars else 2.5
hpad = 0.25
figure_width = nvars*panel_width + (nvars - 1)*hpad
figure_height = panel_height
plt.close(ifig)
fig = plt.figure(ifig, figsize=(figure_width, figure_height), dpi=dpi)
vidx = 0
for v in vars:
print('max(abs({:s})): {:.2e}'.format(v, np.max(np.abs(data[v]))))
if cbars:
ax = plt.Axes(fig, [vidx*(panel_width + hpad)/figure_width, \
(panel_height - panel_width)/panel_height, \
panel_width/figure_width, \
panel_width/panel_height])
else:
ax = plt.Axes(fig, [vidx*(panel_width + hpad)/figure_width, \
0., \
panel_width/figure_width, \
1.])
ax.set_axis_off()
fig.add_axes(ax)
nx, ny, nz = data[v].shape
if axis == 0:
img = np.sum(data[v][(nx//2-slice_width):(nx//2+slice_width), :, :], axis=0)/\
float(2.*slice_width)
elif axis == 1:
img = np.sum(data[v][:, (ny//2-slice_width):(ny//2+slice_width), :], axis=1)/\
float(2.*slice_width)
else:
img = np.sum(data[v][:, :, (nz//2-slice_width):(nz//2+slice_width)], axis=2)/\
float(2.*slice_width)
norm = colors.Normalize(vmin=vlims[v][0], vmax=vlims[v][1], clip=True)
ax.imshow(img, norm=norm, interpolation=interpolation, cmap=cmaps[v])
if v in labels:
for i in range(len(labels[v])):
x, y = (0.5, 0.5)
if label_coords is not None and v in label_coords:
x, y = label_coords[v][i]
color = 'k'
if label_colors is not None and v in label_colors:
color = label_colors[v][i]
ax.text(x, y, labels[v][i], transform=ax.transAxes, color=color)
if cbars:
cax = plt.Axes(fig, [0.05*panel_width/figure_width + vidx*(panel_width + hpad)/figure_width, \
0.5*(panel_height - panel_width)/panel_height, \
0.9*panel_width/figure_width, \
0.5*(panel_height - panel_width)/panel_height])
fig.add_axes(cax)
plt.sca(cax)
cbar = np.transpose(np.repeat(np.linspace(vlims[v][0], vlims[v][1], 256), 16).\
reshape((256, 16)))
extent = [vlims[v][0], vlims[v][1], 0., 1.]
aspect = 0.05*(extent[1] - extent[0])/(extent[3] - extent[2])
cax.imshow(cbar, cmap=cmaps[v], interpolation=interpolation, vmin=vlims[v][0], \
vmax=vlims[v][1], extent=extent, aspect=aspect)
if cbar_labels is not None and v in cbar_labels:
cax.set_xlabel(cbar_labels[v])
cax.get_yaxis().set_visible(False)
cax.get_xaxis().set_tick_params(direction='in')
cax.get_xaxis().set_ticks_position('bottom')
vidx += 1
In [3]:
slice_width = 5
vars = ('Mar', 'Mat', 'rho_fluct')
cmaps = {'Mar':stretch_colormap(cm.viridis, gamma=0.75), \
'Mat':stretch_colormap(cm.viridis, gamma=0.75), \
'rho_fluct':stretch_colormap(cm.BrBG, gamma=0.75, midpoint=0.5)}
cbar_labels = {'Mar':r'Ma$_\mathrm{r}$', \
'Mat':r'Ma$_\perp$', \
'rho_fluct':r'($\rho$ - $\langle\rho\rangle$) / $\langle\rho\rangle$'}
label_coords = {'Mar':((0.025, 0.95), (0.8, 0.95)), \
'Mat':((0.025, 0.95), (0.8, 0.95)), \
'rho_fluct':((0.025, 0.95), (0.8, 0.95))}
label_colors = {'Mar':('w', 'w'), \
'Mat':('w', 'w'), \
'rho_fluct':('k', 'k')}
In [7]:
rid = 'I11'
vlims = {'Mar':(0., 0.2), \
'Mat':(0., 0.2), \
'rho_fluct':(-0.1, 0.1)}
axis = 0
ifig = 1
#for dump in (68, 74, 81, 88):
for dump in (81, ):
t = yp[rid].get('t', fname=dump)[-1]
labels = {'Mar':(rid, '{:.1f} min'.format(t/60.)), \
'Mat':(rid, '{:.1f} min'.format(t/60.)), \
'rho_fluct':(rid, '{:.1f} min'.format(t/60.))}
cbars = True #if dump == 81 else False
plot_panels(moms[rid], dump, slice_width, vars, vlims, cmaps, axis=axis, ifig=ifig, \
cbars=cbars, cbar_labels=cbar_labels, labels=labels, \
label_coords=label_coords, label_colors=label_colors)
plt.savefig('Mar_Mat_rho_fluct_{:s}_{:04d}.pdf'.format(rid, dump), dpi=200)
ifig += 1
In [8]:
rid = 'I5'
vlims = {'Mar':(0., 0.04), \
'Mat':(0., 0.04), \
'rho_fluct':(-0.001, 0.001)}
axis = 0
ifig = 10
#for dump in (130, 260):
for dump in (260, ):
t = yp[rid].get('t', fname=dump)[-1]
labels = {'Mar':(rid, '{:.1f} min'.format(t/60.)), \
'Mat':(rid, '{:.1f} min'.format(t/60.)), \
'rho_fluct':(rid, '{:.1f} min'.format(t/60.))}
cbars = True #if dump == 260 else False
plot_panels(moms[rid], dump, slice_width, vars, vlims, cmaps, axis=axis, ifig=ifig, \
cbars=cbars, cbar_labels=cbar_labels, labels=labels, \
label_coords=label_coords, label_colors=label_colors)
plt.savefig('Mar_Mat_rho_fluct_{:s}_{:04d}.pdf'.format(rid, dump), dpi=150)
ifig += 1
In [ ]: