In [ ]:
%matplotlib inline
import yt
import logging
logging.getLogger('yt').setLevel(logging.ERROR)
import numpy as np
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
import multiprocessing
import sys
sys.path.append('/home/ychen/lib/util')
import util
from tools import setup_cl
#dss = yt.load('/home/ychen/d9/FLASH4/stampede/0529_L45_M10_b1_h1/MHD_Jet_hdf5_plt_cnt*')
In [ ]:
def rescan(dir, printlist=False):
files = util.scan_files(dir, '*hdf5_plt_cnt_[0-9][0-9][0-9][0-9]', \
printlist=printlist, walk=True)
return files
dir = '/home/ychen/d9/FLASH4/stampede/0529_L45_M10_b1_h1/data/'
files = rescan(dir, printlist=True)
In [ ]:
ds = yt.load(files[100].fullpath)
ad = ds.all_data()
zlim = (ds.domain_left_edge[2]/8.0, ds.domain_right_edge[2]/8.0)
print zlim[0].in_units('kpc')
In [ ]:
#sp = ds.sphere([0.0,0.0,0.0], (60, 'kpc'))
prof = yt.create_profile(ad, 'z', ['jet ', 'pressure', 'temperature'], \
logs={'z': False}, weight_field=None, extrema={'z': zlim})
prof_wf = yt.create_profile(ad, 'z', ['jet ', 'pressure', 'temperature'], \
logs={'z': False}, weight_field='cell_mass', extrema={'z': zlim})
In [ ]:
fig, ax1 = plt.subplots()
ax1.plot(prof.x, prof['jet '], c='r', label='pressure')
ax1.plot(prof_wf.x, prof_wf['jet '], c='b', label='pressure')
ax1.semilogy()
ax1.set_ylim(1E-11, 1E6)
In [ ]:
dirs = [
'/home/ychen/data/0only_1110_h0_rerun/',\
'/home/ychen/data/0only_0602_hydro/',\
'/home/ychen/data/0only_0529_h1/',\
'/home/ychen/data/0only_0605_hinf/']
# '/home/ychen/data/0only_1212_h0_10Myr_rerun/',\
# '/home/ychen/data/0only_0413_hydro_10Myr/',\
# '/home/ychen/data/0only_1022_h1_10Myr/',\
# '/home/ychen/data/0only_0204_hinf_10Myr/']
colors, labels = setup_cl(dirs)
tt = np.linspace(0, 5, 100)
for dir in dirs:
data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
label = labels[dir]
if '10Myr' in dir:
plt.fill_between(data[:,1], data[:,2], -data[:,3],\
facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.8, label=label)
#plt.plot(data[:,1], data[:,2], c=colors[dir])
else:
plt.fill_between(data[:,1], data[:,2], -data[:,3],\
facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.8, label=label)
#plt.fill_between(data[:,1], data[:,4], -data[:,5],\
# facecolor=colors[dir], lw=1, linestyle=':', edgecolor=colors[dir], alpha=0.3)
plt.plot(tt, 10*tt, ':', c='orange', label='10 kpc / Myr')
plt.plot(tt, 8*tt, '--', c='cyan', label='8 kpc / Myr')
plt.legend(loc=4)
handles, labels = plt.axes().get_legend_handles_labels()
# reverse the order
plt.legend(handles[::-1], labels[::-1], loc=4)
plt.grid(ls='--', alpha=0.5)
plt.ylabel('z (kpc)')
plt.xlabel('t (Myr)')
plt.xlim(0,20)
plt.ylim(0,50)
plt.savefig('LobeSize_20Myr.pdf')
#plt.semilogx()
#plt.semilogy()
In [ ]:
dirs = [
# '/home/ychen/data/0only_1110_h0_rerun/',\
# '/home/ychen/data/0only_0602_hydro/',\
# '/home/ychen/data/0only_0529_h1/',\
# '/home/ychen/data/0only_0605_hinf/',\
'/home/ychen/data/0only_1212_h0_10Myr_rerun/',\
'/home/ychen/data/0only_0413_hydro_10Myr/',\
'/home/ychen/data/0only_1022_h1_10Myr/',\
'/home/ychen/data/0only_0204_hinf_10Myr/']
colors, labels = setup_cl(dirs)
tt = np.linspace(0, 5, 100)
for dir in dirs:
data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
label = labels[dir]
if '10Myr' in dir:
plt.fill_between(data[:,1], data[:,2], -data[:,3],\
facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.8, label=label)
#plt.plot(data[:,1], data[:,2], c=colors[dir])
else:
plt.fill_between(data[:,1], data[:,2], -data[:,3],\
facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.8, label=label)
#plt.fill_between(data[:,1], data[:,4], -data[:,5],\
# facecolor=colors[dir], lw=1, linestyle=':', edgecolor=colors[dir], alpha=0.3)
#plt.plot(tt, 10*tt, ':', c='orange', label='10 kpc / Myr')
#plt.plot(tt, 8*tt, '--', c='cyan', label='8 kpc / Myr')
plt.legend(loc=4)
handles, labels = plt.axes().get_legend_handles_labels()
# reverse the order
plt.legend(handles[::-1], labels[::-1], loc=4)
plt.ylabel('z (kpc)')
plt.xlabel('t (Myr)')
plt.xlim(0.2,100)
plt.ylim(2,80)
plt.semilogx()
plt.semilogy()
xticks = [0.3, 1,3,10,30,100]
plt.xticks(xticks, xticks)
yticks = [3,10,30,80]
plt.yticks(yticks, yticks)
plt.grid(ls='--', alpha=0.5)
plt.savefig('LobeSize_100Myr.pdf')
In [ ]:
data[:3]
In [ ]:
dirs = [
'/home/ychen/data/0only_0605_h0/',\
'/home/ychen/data/0only_1110_h0_rerun/']
# Set up colors and label names
colors = {}
labels = {}
for dirname in dirs:
if '0605_h0' in dirname:
colors[dirname] = 'b'
labels[dirname] = 'poloidal'
elif '1110_h0_rerun' in dirname:
colors[dirname] = 'cyan'
labels[dirname] = 'poloidal (rerun)'
fig = plt.figure(figsize=(8,6))
for dir in dirs:
data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
label = labels[dir]
tt = data[:,1]
zEdge1 = data[:,2]
zEdge2 = -data[:,3]
rEdge1 = data[:,6]
rEdge2 = data[:,7]
plt.fill_between(tt, zEdge1, zEdge2,\
facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.7, label=label+' zEdge')
plt.fill_between(tt, rEdge1, rEdge2,\
facecolor=colors[dir], lw=1, linestyle=':', edgecolor=colors[dir], alpha=0.3, label=label+' rEdge')
plt.legend(loc=2)
plt.ylabel('Edge (kpc)')
plt.xlabel('t (Myr)')
plt.xlim(0,20)
plt.ylim(0,50)
#plt.semilogx()
In [ ]:
dirs = [
'/home/ychen/data/0only_0602_hydro/',\
'/home/ychen/data/0only_0529_h1/',\
'/home/ychen/data/0only_0605_hinf/',\
'/home/ychen/data/0only_1106_M3_h1/',\
'/home/ychen/data/0only_1204_M24_b01/',\
'/home/ychen/data/0only_1110_h0_rerun/']
fig = plt.figure(figsize=(8,6))
for dir in dirs:
data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
#label = dir[28:-1]
label = labels[dir]
plt.fill_between(data[:,1], data[:,2], -data[:,3],\
facecolor=colors[dir], lw=1, edgecolor=colors[dir], alpha=0.8, label=label)
#plt.fill_between(data[:,1], data[:,4], -data[:,5],\
# facecolor=colors[dir], lw=1, linestyle=':', edgecolor=colors[dir], alpha=0.3)
plt.legend(loc=2)
plt.ylabel('zEdge (kpc)')
plt.xlabel('t (Myr)')
plt.savefig('LobsSize_and_LobeCenter.pdf')
#plt.axvline(10, ls=':', c='grey')
In [ ]:
dirs = [
'/home/ychen/data/0only_0602_hydro/',\
'/home/ychen/data/0only_0529_h1/',\
'/home/ychen/data/0only_0605_hinf/',\
'/home/ychen/data/0only_1106_M3_h1/',\
'/home/ychen/data/0only_1204_M24_b01/',\
'/home/ychen/data/0only_1110_h0_rerun/']
fig = plt.figure(figsize=(8,6))
for dir in dirs:
data = np.loadtxt(dir+'GridAnalysis_LobeSize.txt')
#label = dir[28:-1]
label = labels[dir]
tt = data[:,1]
zEdge1 = data[:,2]
zEdge2 = -data[:,3]
zEdge = (zEdge1+zEdge2)/2.0
rEdge1 = data[:,6]
rEdge2 = data[:,7]
rEdge = (rEdge1+rEdge2)/2.0
zrRatio = zEdge/rEdge
plt.plot(tt, zrRatio, lw=1, color=colors[dir], label=label)
plt.legend(loc=1)
plt.ylabel('zEdge/rEdge')
plt.xlabel('t (Myr)')
plt.ylim(2, 5)
#plt.savefig('LobsSize_and_LobeCenter.pdf')
#plt.axvline(10, ls=':', c='grey')
In [ ]:
#upfilter = ad['z'] > 0.0
kpc = yt.units.kpc.in_units('cm')
leftedge = [-30*kpc, -30*kpc, 0.0]
rightedge = [ 30*kpc, 30*kpc, 60*kpc]
upbox = ds.box(leftedge, rightedge)
In [ ]:
In [ ]:
np.sum(upbox['jet ']*upbox['cell_mass']*upbox['z'])/np.sum(upbox['jet ']*upbox['cell_mass'])
In [ ]:
fig, ax1 = plt.subplots()
ax1.plot(prof_p.x, prof_p['pressure'], c='b', label='pressure')
plt.legend(loc=2)
ax2 = ax1.twinx()
ax2.plot(prof_t.x, prof_t['temperature'], c='r', label='temperature')
plt.legend(loc=1)
#plot = yt.ProfilePlot.from_profiles([prof_p, prof_t], labels=['pressure', 'temperature'])
In [ ]:
plt.plot(prof[0].x[1:], (prof[0]['temperature'][1:]-prof[0]['temperature'][0:-1])/)
In [ ]:
In [ ]: