Mach number and density fluctuations in runs I11 and I5


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)


Reading attributes from file  YProfile-01-0266.bobaaa
There are 267 YProfile files in the /data/ASDR/PPMstar/C-ingestion/I5/YProfile/ directory.
Ndump values range from 1 to 267
Time values range from 8.87622 to 2369.95
Reading attributes from file  YProfile-01-0090.bobaaa
There are 91 YProfile files in the /data/ASDR/PPMstar/C-ingestion/I11/YProfile/ directory.
Ndump values range from 1 to 91
Time values range from 8.87622 to 807.736

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


max(abs(Mar)): 2.67e-01
max(abs(Mat)): 2.87e-01
max(abs(rho_fluct)): 1.48e-01

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


max(abs(Mar)): 4.45e-02
max(abs(Mat)): 5.39e-02
max(abs(rho_fluct)): 4.04e-02

In [ ]: