In [ ]:
%matplotlib inline
import numpy as np
import util
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
matplotlib.rcParams['animation.html'] = 'html5'
import matplotlib.pyplot as plt
import h5py
import yt
from tools import calcDen0_2015
from particle_filters import deltaden0
me = yt.utilities.physical_constants.mass_electron.v
c = yt.utilities.physical_constants.speed_of_light.v
e = yt.utilities.physical_constants.elementary_charge.v
Myr= yt.units.Myr.in_units('s').v
kpc= yt.units.kpc.in_units('cm').v
In [ ]:
_fields_list = [
# 'tag', 'tadd', 'den0',
'posx', 'posy', 'posz',
'velx', 'vely', 'velz',
'dens', 'gamc', 'tau',
'magx', 'magy', 'magz',
'shok'
]
#dir = '/home/ychen/d9/FLASH4/2015_production_runs/0529_L45_M10_b1_h1/'
#plotfiles = util.scan_files(dir, 'MHD_Jet_hdf5_plt_cnt_00[0-9][0-9]', printlist=False)
#partfiles = util.scan_files(dir, 'MHD_Jet_hdf5_part_????_updated', printlist=False)
def find_part_ind(h5file, (tadd, tag)):
"""
Find the index of the particle matching tadd and tag.
"""
tp = h5file['tracer particles']
colname = [item[0].strip() for item in h5file['particle names']]
itadd = colname.index('tadd')
itag = colname.index('tag')
ind_tadd = np.isclose(tp[:,itadd]-tadd, 0.0)
ind_tag = np.in1d(tp[:,itag], tag)
ind = np.where(np.logical_and(ind_tadd, ind_tag))[0]
#print np.where(np.in1d(tp[:,findices['tadd']], tadd))
if len(ind) == 0:
return None
elif len(ind) == 1:
return ind
else:
raise IndexError
def find_part_indices(h5file, tags, tadds):
"""
Find the indices of the particles matching tadds and tags. Return None if not found.
"""
tp = h5file['tracer particles']
colname = [item[0].strip() for item in h5file['particle names']]
itadd = colname.index('tadd')
itag = colname.index('tag')
masks_tadd = np.in1d(tp[:,itadd], tadds)
masks_tag = np.in1d(tp[:,itag], tags)
indices_temp = np.where(np.logical_and(masks_tadd, masks_tag))[0]
assert len(indices_temp) <= len(tags)
indices = [None]*len(tags)
for i, (tag, tadd) in enumerate(zip(tags, tadds)):
for ind in indices_temp:
if tag == tp[ind,itag] and np.isclose(tadd, tp[ind,itadd]):
indices[i] = ind
break
return indices
class Particle():
def __init__(self, (tag, tadd), grid_fields=None):
self.tadd = tadd
self.tag = tag
self.den0 = -1
self.time = []
#nfiles = len(partfiles)
self.grid_fields = grid_fields
if grid_fields:
for grid_field in grid_fields:
setattr(self, grid_field, [])
for field in _fields_list:
setattr(self, field, [])
def read_from_h5file(self, ind, h5file, ds=None):
# If the supplied particle is found in this particle file
if ind:
self.time.append( h5file['real scalars'].value[0][1]/Myr )
tp = h5file['tracer particles']
# Assign den0 for the first time and make sure den0 are the same
if self.den0 < 0:
self.den0 = tp[ind, colname.index('den0')]
else:
assert self.den0 == tp[ind, colname.index('den0')]
for field in _fields_list:
getattr(self, field).append( tp[ind,findices[field]] )
if self.grid_fields and ds:
for grid_field in self.grid_fields:
pos = [self.posx[-1], self.posy[-1], self.posz[-1]]
#print ds, grid_field, pos
getattr(self, grid_field).append(ds.find_field_values_at_point(grid_field, pos))
class Particles(Particle):
def __init__(self, tags, tadds, grid_fields=None):
assert len(tadds) == len(tags)
self.tags = tags
self.tadds = tadds
self.nparticles = len(tadds)
self.grid_fields = grid_fields
# Initialize the particle
self.data = [Particle((tag,tadd), grid_fields=grid_fields) for tag, tadd in zip(tags, tadds)]
def __getitem__(self, key):
return self.data[key]
def read_from_partfile(self, pf):
print pf
h5file = h5py.File(pf.fullpath, 'r')
colname = [item[0].strip() for item in h5file['particle names']]
indices = find_part_indices(h5file, self.tags, self.tadds)
for i, ind in enumerate(indices):
ds = yt.load(pf.fullpath.replace('part', 'plt_cnt')) if self.grid_fields else None
self.data[i].read_from_h5file(ind, h5file, ds=ds)
def read_from_partfiles(self, partfiles):
for pf in partfiles:
self.read_from_partfile(pf)
# Convert the lists to numpy arrays
for part in self.data:
for field in _fields_list:
setattr(part, field, np.array(getattr(part, field)))
if self.grid_fields:
for grid_field in self.grid_fields:
setattr(part, grid_field, np.array(getattr(part, grid_field)))
In [ ]:
dir = '/d/d9/ychen/FLASH4/2015_production_runs/0529_L45_M10_b1_h1/'
partfiles = util.scan_files(dir, 'MHD_Jet_hdf5_part_??00', printlist=True)
h5file = h5py.File(partfiles[6].fullpath, 'r')
colname = [item[0].strip() for item in h5file['particle names']]
findices = {f: colname.index(f) for f in _fields_list}
# Arbitrarily pick a particle
import random
tp = h5py.File(partfiles[6].fullpath, 'r')['tracer particles']
nparticles = 10
rints = sorted([random.randint(0, tp.shape[0]) for i in range(nparticles)])
tags = tp[rints,15]
tadds = tp[rints,8]
shoks = tp[rints,7]
for (tag, tadd, shok, ind) in zip(tags, tadds, shoks, rints):
print tag, tadd, shok, ind
#tadd, tag = neg_dens_list[1600]
#tag, tadd = 127225.0, 218857882571403.5
#tag, tadd = 14101.0, 14079945645129.127
#part = Particle((tag,tadd), partfiles, grid_fields=['pressure'])
In [ ]:
parts = Particles(tags, tadds, grid_fields=['pressure'])
parts.read_from_partfiles(partfiles)
In [ ]:
for part in parts:
print part.tag, part.tadd, len(part.posx)
In [ ]:
tp = h5py.File(partfiles[50].fullpath, 'r')['tracer particles']
indices_tadd = np.in1d(tp[:,8], tadds)
indices_tag = np.in1d(tp[:,15], tags)
indices = find_part_indices(h5py.File(partfiles[90].fullpath, 'r'), tags, tadds)
for (tag, tadd, ind) in zip(tags, tadds, indices):
print tag, tadd, ind
In [ ]:
import pickle
parts = pickle.load(open("Particles.pickle", "r"))
line = np.linspace(-13,-7)
for part in parts[:1]:
size = 9-part.grid_level
fig = plt.figure()
ax = fig.add_subplot(111)
E_B = (part.magx**2 + part.magy**2 + part.magz**2)/2
Br = np.sqrt(part.magx**2 + part.magy**2)
v = np.log10(np.sqrt(part.velx**2 + part.vely**2 + part.velz**2)) - 5
r = np.sqrt(part.posx**2 + part.posy**2)/kpc
sc = ax.scatter(r, np.abs(part.posz)/kpc, lw=0, s=2, vmin=2.5, vmax=4.5, c=np.log10(part.entropy))
ax.set_xlim(0,40)
ax.set_xlabel('r (kpc)')
ax.set_ylim(0,40)
ax.set_ylabel('|z| (kpc)')
ind = -1
fieldqv = ax.quiver(r[ind], np.abs(part.posz[ind])/kpc, Br[ind], np.abs(part.magz[ind]), \
headwidth=0, headlength=0, headaxislength=0, width=0.002, pivot='middle', scale=1E-4)
cb = plt.colorbar(sc, pad=0.1)
cb.set_label('log entropy | log v (km/s)')
txt = ax.text(30, 2, '%.2f Myr' % part.time[-1])
ax2 = plt.twinx().twiny()
ax2.plot(line, line, c='gray')
ax2.plot(line, line-1, c='darkgray')
ax2.plot(line, line-2, c='silver')
ax2.plot(line, line-3, c='lightgray')
plt.plot(np.log10(part.pressure), np.log10(E_B), lw=1, alpha=0.1)
#plt.scatter(np.log10(part.pressure), np.log10(E_B), lw=0, c=part.grid_level, vmin=5, vmax=8, s=5)
sc2 = plt.scatter(np.log10(part.pressure), np.log10(E_B), c=v, vmin=2.5, vmax=4.5, lw=0, s=size*3)
txt = '%.2f Myr, %i' % (part.tadd/Myr, part.tag)
if part.shok[0] > 0.5:
txt += ' shock'
plt.text(-12.5,-7.5, txt )
ax2.set_xlim(-13, -7)
ax2.set_xlabel('log(P_gas)')
ax2.set_ylim(-13, -7)
ax2.set_ylabel('log(P_B)')
fname = 'particle_%.2fMyr_%i.png' % (part.tadd/Myr, part.tag)
#plt.savefig('/home/ychen/www/MHD_Jet/201705/P_PB_v/' + fname)
#plt.show()
In [ ]:
import pickle
fname = '/d/d11/ychen/2015_production_runs/1022_L45_M10_b1_h1_10Myr/Particles_tracing.pickle'
parts = pickle.load(open(fname, "rb"))
line = np.linspace(-13,-7)
for part in parts[:1]:
size = 9-part.grid_level
plt.figure()
E_B = (part.magx**2 + part.magy**2 + part.magz**2)/2
beta = np.array([p/pb for (p,pb) in zip(part.pressure,E_B)])
v = np.log10(np.sqrt(part.velx**2 + part.vely**2 + part.velz**2))
plt.scatter(part.grid_level, np.log10(beta), c='r', lw=0, s=size*3)
#plt.scatter(part.time, np.log10(E_B), c='b', lw=0, s=size*3)
txt = '%.2f Myr, %i' % (part.tadd/Myr, part.tag)
#if part.shok[0] > 0.5:
# txt += ' shock'
#plt.text(-12.5,-7.5, txt )
#cb=plt.colorbar()
#cb.set_label(u'log $v$')
#plt.ylim(-13, -7)
#plt.ylabel('log(P)')
#plt.ylim(0,100)
#plt.twinx()
#plt.plot(part.time, part.grid_level, c='gray')
#plt.xlim(0, 25)
#plt.xlabel('time (Myr)')
#fname = 'particle_%.2fMyr_%i.png' % (part.tadd/Myr, part.tag)
#plt.savefig('/home/ychen/www/MHD_Jet/201705/P_PB_v/' + fname)
plt.show()
In [ ]:
from tempfile import NamedTemporaryFile
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
import numpy as np
import base64
IMG_TAG = """<img src="data:image/gif;base64,{0}" alt="some_text">"""
def anim_to_gif(anim):
data="0"
with NamedTemporaryFile(suffix='.gif') as f:
anim.save(f.name, writer='imagemagick', fps=10);
data = open(f.name, "rb").read()
data = base64.b64encode(data)
return IMG_TAG.format(data.decode())
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_gif(anim))
fig = plt.figure()
fig.set_dpi(100)
fig.set_size_inches(7, 6.5)
ax = plt.axes(xlim=(0, 10), ylim=(0, 10))
patch = plt.Circle((5, -5), 0.75, fc='y')
patch2 = plt.Circle((5, -5), 0.1, fc='b')
c= np.array([5,5])
import pickle
with open("Particles.pickle", "rb") as f:
parts = pickle.load(f, encoding='latin1')
eqline = np.linspace(-13,-7)
for part in parts:
size = 9-part.grid_level
fig = plt.figure()
ax = fig.add_subplot(111)
E_B = (part.magx**2 + part.magy**2 + part.magz**2)/2
Br = np.sqrt(part.magx**2 + part.magy**2)
v = np.log10(np.sqrt(part.velx**2 + part.vely**2 + part.velz**2)) - 5
r = np.sqrt(part.posx**2 + part.posy**2)
entr = np.log10(part.entropy)
sc = ax.scatter(r/kpc, np.abs(part.posz)/kpc, lw=0, s=1, c=entr, vmin=2.5, vmax=4.5)
ax.set_xlim(0,40)
ax.set_xlabel('r (kpc)')
ax.set_ylim(0,40)
ax.set_ylabel('|z| (kpc)')
cb = plt.colorbar(sc, pad=0.1)
cb.set_label(u'left: log entropy | right: log $v$ (km/s)')
time_txt = plt.text(25, 2, '%.2f Myr' % part.time[0])
fieldqv = ax.quiver(r[0]/kpc, np.abs(part.posz[0])/kpc, Br[0], np.abs(part.magz[0]), \
headwidth=0, headlength=0, headaxislength=0, width=0.002, pivot='middle', scale=1E-4)
ax2 = plt.twinx().twiny()
ax2.plot(eqline, eqline, c='gray')
ax2.plot(eqline, eqline-1, c='darkgray')
ax2.plot(eqline, eqline-2, c='silver')
ax2.plot(eqline, eqline-3, c='lightgray')
line, = ax2.plot(np.log10(part.pressure), np.log10(E_B), lw=1, alpha=0.1)
#plt.scatter(np.log10(part.pressure), np.log10(E_B), lw=0, c=part.grid_level, vmin=5, vmax=8, s=5)
sc2 = ax2.scatter(np.log10(part.pressure), np.log10(E_B), c=v, vmin=2.5, vmax=4.5, lw=0, s=2)
txt = '%.2f Myr, %i' % (part.tadd/Myr, part.tag)
if part.shok[0] > 0.5:
txt += ' shock'
ax2.text(-12.5,-7.5, txt )
ax2.set_xlim(-13, -7)
ax2.set_xlabel('log(P_gas)')
ax2.set_ylim(-13, -7)
ax2.set_ylabel('log(P_B)')
#plt.savefig('/home/ychen/www/MHD_Jet/201705/P_PB_v/' + fname)
#plt.show()
def animate(i):
line.set_data(np.log10(part.pressure[:i]), np.log10(E_B[:i]))
#plt.scatter(np.log10(part.pressure), np.log10(E_B), lw=0, c=part.grid_level, vmin=5, vmax=8, s=5)
sc.set_offsets((r[:i]/kpc, np.abs(part.posz[:i])/kpc))
sc.set_array(entr[:i])
sc2.set_offsets((np.log10(part.pressure[:i]), np.log10(E_B[:i])))
sc2.set_array(v[:i])
time_txt.set_text('%5.2f Myr, res:%i' % (part.time[i-1], part.grid_level[i-1]))
fieldqv.set_offsets((r[i-1]/kpc, np.abs(part.posz[i-1])/kpc))
fieldqv.set_UVC(Br[i-1], np.abs(part.magz[i-1]))
#sc.set_sizes(size[:i]*2)
return line, sc, sc2, time_txt, fieldqv
anim = animation.FuncAnimation(fig, animate, frames=len(part.magx), interval=1)
#fname = 'particle_%.2fMyr_%i.mp4' % (part.tadd/Myr, part.tag)
#anim.save('/home/ychen/www/MHD_Jet/201705/P_PB_v/' + fname, writer='ffmpeg', fps=10, bitrate=8E6)
display_animation(anim)
#anim.to_html5_video()
In [ ]:
v_thres = 0.5*3E9
def plot_particle_tracing(part, timerange='all'):
if timerange == 'all':
part.time -= part.time[0]
part.time[0] += 1E-2
fig = plt.figure(figsize=(10,6))
#################### Density ################
#rho_norm = 10**int(np.log10(max(part.dens)))
#rho_norm = max(part.dens)/10
rho_norm = 1E-26
plt.fill_between(part.time, part.dens/rho_norm, lw=0, color='skyblue', alpha=0.7,
label=u'density/(%.0eg/cm$^3$)' % rho_norm)
plt.fill_between(part.time, part.den1/rho_norm, lw=0, color='skyblue', alpha=0.5)
#plt.plot(part.time, part.tau, '^-', c='r')
#dtau = np.array(part.tau[1:])-np.array(part.tau[:-1])
#plt.plot(part.time, dtau, '-', c='r')
#################### Velocity ################
v_mag = np.sqrt(part.velx**2+part.vely**2+part.velz**2)
v_norm = 3E8
#plt.plot(part.time, v_mag/v_norm, '--', c='r', alpha=0.5, label='v/%.0e cm/s' % v_norm)
plt.fill_between(part.time, v_mag/v_norm, color='r', lw=0, alpha=0.3, label=u'$|v| & v_z $/(%.0ecm/s)' % v_norm)
plt.fill_between(part.time, np.abs(part.velz)/v_norm, color='r', lw=0, alpha=0.2)
arg = np.where(v_mag<v_thres)[0][0]
plt.axvline(part.time[arg], color='r', lw=1, alpha=0.3)
#plt.hlines(v_thres/v_norm, part.time[0], part.time[-1], color='r', lw=1, alpha=0.3)
#plt.twiny()
#################### Magnetic Field ################
B_mag = np.sqrt(part.magx**2+part.magy**2+part.magz**2)
#B_norm = max(B_mag)/10
B_norm = 1E-5
plt.plot(part.time, B_mag/B_norm, '-', c='b', lw=0.8, alpha=0.5, label=u'B/(%.0f$\mu$G)' % (B_norm*1E6))
#################### TAU ################
tau_norm = max(part.tau)/10
plt.plot(part.time, part.tau/tau_norm, ':', c='g', label='tau/%.0e' % (tau_norm), alpha=0.3)
plt.plot(part.time, part.dtau/tau_norm, '--', c='g', label='dtau/%.0e' % (tau_norm))
#plt.hlines(0.00056966283093648636/tau_norm, part.time[0], part.time[-1], color='g', lw=0.5)
#ind = np.argwhere(np.isclose(part.tau, 0.00056966283093648636))
#plt.vlines(part.time[ind], 0, 12, color='g', lw=0.5)
#################### Cutoff Gamma ################
plt.plot(part.time, np.log10(part.gamc), ':', c='c', label=u'log $\\gamma_{c}$', alpha=0.3)
gamc_dtau = (part.dens/part.den1)**(1./3.)/part.dtau
plt.plot(part.time, np.log10(gamc_dtau), '-', c='c', label=u'log $\\gamma_{c}$ (dtau)')
#################### Cutoff frequency ################
nuc = 3.0*part.gamc**2*e*B_mag/(4.0*np.pi*me*c)
nuc_dtau = 3.0*gamc_dtau**2*e*B_mag/(4.0*np.pi*me*c)
plt.plot(part.time, np.log10(nuc), ':', c='gold', label=u'log $\\nu_{c}$/Hz', alpha=0.3)
plt.plot(part.time, np.log10(nuc_dtau), '-', c='gold', label=u'log $\\nu_{c}$/Hz (dtau)')
#################### Cylindrical Radius ################
#r_nozzle = 7.5E20
#rr = np.sqrt(part.posx**2 + part.posy**2)/r_nozzle
#plt.plot(part.time, rr, '-', c='indigo', label=u'$r_{cylindrical}$/(%ipc)' % (r_nozzle/kpc*1000))
plt.plot(part.time, [14]*len(part.time), '|', c='k')
plt.ylim(0,14)
plt.grid(axis='y', alpha=0.5)
if timerange == 'all':
plt.xlim(part.time[0], part.time[-1])
plt.xlabel('time (Myr) - t0')
plt.semilogx()
else:
plt.xlim(part.time[0], part.time[0]+timerange)
plt.xlabel('time (Myr)')
plt.legend(loc=1, ncol=2)
plt.text(0.05, 0.92,
't0=%.2f Myr, id=%i' % (part.tadd/Myr, part.tag),
transform=plt.axes().transAxes)
plt.text(0.05, 0.87,
r'$\rho_1$=%.2f*%.0e g/cm$^3$' % (part.den1/rho_norm, rho_norm),
transform=plt.axes().transAxes)
coreden0 = calcDen0_2015(part.tadd)
plt.plot(part.time[0], coreden0/rho_norm, '+', c='skyblue', ms=5)
if part.shok[0] == 1:
plt.text(0.05, 0.82, 'shok', transform=plt.axes().transAxes)
elif part.den0 < coreden0 - deltaden0 or part.den0 > coreden0 + deltaden0:
plt.text(0.05, 0.82, r'$\rho_0$=%.3e, $\rho_{0,core}=$%.3e' % (part.den0, coreden0), transform=plt.axes().transAxes)
plt.text(0.05, 0.77, 'shelth', transform=plt.axes().transAxes)
else:
plt.text(0.05, 0.82, r'$\rho_0$=%.3e, $\rho_{0,core}=$%.3e' % (part.den0, coreden0), transform=plt.axes().transAxes)
if timerange == 'all':
figfname = '../figures/particle_tracing/alltime/particle_tracing_%.2f_%i_%.2f.png'
else:
figfname = '../figures/particle_tracing/%iMyr/particle_tracing_%%.2f_%%i_%%.2f.png' % timerange
plt.tight_layout()
#plt.savefig(figfname % (part.tadd/Myr, part.tag, part.den1/rho_norm))
#plt.close(fig)
plt.show()
In [ ]:
fname = '/d/d11/ychen/2015_production_runs/1022_L45_M10_b1_h1_10Myr/Particles_tracing.pickle'
parts = pickle.load(open(fname, "rb"))
for i, part in enumerate(parts[-1:]):
plot_particle_tracing(part, timerange=2)
plot_particle_tracing(part)
In [ ]:
neg_dens_list = []
def get_field_from_h5file(h5file, field):
tp = h5file['tracer particles']
colname = [item[0].strip() for item in h5file['particle names']]
ifield = colname.index(field)
return tp[:,ifield]
for pfile in partfiles:
h5file = h5py.File(pfile.fullpath, 'r')
denss = get_field_from_h5file(h5file, 'dens')
tadds = get_field_from_h5file(h5file, 'tadd')
tags = get_field_from_h5file(h5file, 'tag')
mask = np.where(denss < 0)
#print np.count_nonzero(mask)
neg_dens_list.extend( zip(tadds[mask], tags[mask]) )
plt.scatter(denss[mask], tadds[mask])
plt.xlim(-1E-26,0)
In [ ]:
print len(neg_dens_list), len(set(neg_dens_list))
In [ ]:
p = yt.SlicePlot(ds, "x", "pressure", center='c')
p.zoom(1)
p.annotate_particles(slab_width, col='k', marker='o')
part_tag = 7.0
part_ind = abs(ad['particle_tag']-part_tag).argmin()
textfmt = r'$\gamma_c = %11.3e$'+'\n'+r'$\rho = %11.3e$'+'\n'+'$|B| = %11.3e$'
B = sqrt(ad['particle_magx'][part_ind]**2 + ad['particle_magy'][part_ind]**2 + ad['particle_magz'][part_ind]**2)*sqrt(4.0*np.pi)
text = textfmt % (ad['particle_gamc'][part_ind], ad['particle_dens'][part_ind], B)
text_args = {'color': 'w', 'size':16, 'ha':'right'}
marker_args = {'c':'cyan', 's':20}
pos = (ad['particle_position_y'][part_ind], ad['particle_position_z'][part_ind])
p.annotate_marker(pos, marker='o', plot_args=marker_args)
p.annotate_text((0.95,0.89), text, text_args=text_args)
p.show()
In [ ]: