In [ ]:
%matplotlib inline
import h5py
import numpy as np
import matplotlib.pyplot as plt
dir = '/home/ychen/d9/FLASH4/2015_production_runs/0529_L45_M10_b1_h1/'
f = h5py.File(dir+'MHD_Jet_hdf5_part_0600', 'r')
In [ ]:
f.values()
val = f.values()[-2]
In [ ]:
particles = []
#colname = ['den0','dens','gamc','magx','magy','magz','tadd','tau','blk','posx','posy','posz','proc','tag','velx','vely','velz']
colname = [item[0].strip() for item in f.values()[5]]
print colname
for value_particle in val.value:
particle = dict(zip(colname, value_particle))
particles.append(particle)
particles.sort(key=lambda p:p['tag'])
In [ ]:
fields = ['tag', 'proc', 'tadd', 'den0', 'dens', 'gamc', 'tau', 'posz']
print ('%11s'*len(fields)) % tuple(fields)
for par in particles[-7600:-7590]:
B2 = par['magx']*par['magx'] + par['magy']*par['magy'] + par['magz']*par['magz']
strfmt = '%11.3e'*len(fields)
field_values = [par[field] for field in fields]
print strfmt % tuple(field_values)
In [ ]:
f["particle names"][:]
In [ ]:
import h5py
dir = '/home/ychen/d9/FLASH4/2016_test/1028_metal_test/'
f = h5py.File(dir+'MHD_Jet_hdf5_part_0001', 'r')
In [ ]:
for key, values in f.items():
print(key, values)
In [ ]:
print(f['particle names'].value)
print(f['tracer particles'][:100,10])
In [ ]:
ds = f.values()[6]
id = 9193
print 'coordinates = ', f.get('coordinates').value[id]
print 'refine level = ', f.get('refine level').value[id]
print 'which child = ', f.get('which child').value[id]
print f.get('gid').value[id]
print f.get('gsurr_blks').value[id][:,:,:,0]
In [ ]:
tags = np.array([p['tag'] for p in particles])
tadds = np.array([p['tadd'] for p in particles])
denss = np.array([p['dens'] for p in particles])
den0s = np.array([p['den0'] for p in particles])
newtags = np.array([p['tadd']+p['tag'] for p in particles])
In [ ]:
mask = np.where(denss < 0)
print np.count_nonzero(mask)
plt.scatter(denss[mask], tadds[mask])
plt.xlim(-1E-26,0)
In [ ]:
print len(tags)
print len(set(tags))
print len(set(newtags))
In [ ]:
plt.figure(figsize=(6,4))
plt.scatter(tadds, tags, lw=0, s=1)
In [ ]:
import util
import pickle
import h5py
import numpy as np
dir = '/home/ychen/d9/FLASH4/2015_production_runs/0529_L45_M10_b1_h1/'
def rescan(printlist=False):
files = util.scan_files(dir, '*hdf5_part_[0-9][0-9][0-9][0-9]', printlist=printlist)
return files
files = rescan(True)
In [ ]:
# Velocity threshold
v_thres = 0.01*3E10
# Particles that leave the jet
particles_leave = {}
for f in files[:]:
h5f = h5py.File(f.fullpath, 'r')
colname = [item[0].strip() for item in h5f['particle names']]
tag = colname.index('tag')
tadd = colname.index('tadd')
dens = colname.index('dens')
tau = colname.index('tau')
shok = colname.index('shok')
velx = colname.index('velx')
vely = colname.index('vely')
velz = colname.index('velz')
# Fields to be written in the output file
fields = [tadd, tag, dens, tau]
val = h5f['tracer particles']
# Go through the list of particles
for part in val.value:
vx = part[velx]
vy = part[vely]
vz = part[velz]
vel_magnitude = np.sqrt(vx*vx+vy*vy+vz*vz)
#if part[shok] > 0.01:
# if part[shok] < 1.0:
# print par[shok]
# next
# Record new particles that just leave the jet
if (part[tag], part[tadd]) not in particles_leave and vel_magnitude < v_thres:
particles_leave[(part[tag], part[tadd])] = [part[field] for field in fields]
# Collect particle values to a list of dictionaries
#particles = []
#for value_particle in val.value:
# particle = dict(zip(colname, value_particle))
# particles.append(particle)
# Go through the list of particles
#for part in particles:
# newtag = part['tadd']+part['tag']
# if part['shok'] > 0.01:
# if part['shok'] < 1.0:
# print part['shok']
# next
# Record new particles that just leave the jet
# if newtag not in particles_leave and part['velz'] < v_thres:
# particles_leave[newtag] = [part[field] for field in fields]
#print f, '%6i : %6i' % (len(particles), len(particles_leave))
#particles_leave = np.array([[key[0] + key[1]] + value for key, value in particles_leave.iteritems()])
pickle.dump(particles_leave, open('%s_particles_leave.pickle' % dir.split('/')[-2], 'wb'))
In [ ]:
read = pickle.load(open('%s_particles_leave.pickle' % dir.split('/')[-2], 'r'))
print read.shape
In [ ]:
h5f = h5py.File(files[64].fullpath, 'r')
val = h5f['tracer particles']
vx = val.value[:,velx]
vy = val.value[:,vely]
vz = val.value[:,velz]
vel_magnitude = np.sqrt(vx*vx+vy*vy+vz*vz)
print val.value[vel_magnitude < v_thres].shape
read = pickle.load(open('%s_particles_leave.pickle' % dir.split('/')[-2], 'r'))
print read.shape
In [ ]:
h5f = h5py.File(files[64].fullpath, 'r')
val = h5f['tracer particles']
print val.shape
In [ ]:
for f in files[::10]:
print f
h5f = h5py.File(f.fullpath, 'r')
val = h5f['tracer particles']
newtagss = val.value[:,tag]+val.value[:,tadd]
ind = np.where(np.in1d(newtagss, read[:,0], invert=True))[0]
if ind.shape[0]>0:
print ind.shape
np.where(val.value[ind,velz] > v_thres)
In [ ]:
ind = np.where(np.in1d(read[:,2], 12212.0))[0]
read[ind,1]
In [ ]:
f = files[65]
h5f = h5py.File(f.fullpath, 'r')
colname = [item[0].strip() for item in h5f['particle names']]
val = h5f['tracer particles']
In [ ]:
val.value[np.where(np.in1d(val.value[:,tag], 31947.0))[0], tadd]
In [ ]:
ind = np.where(np.in1d(read[:,2], 1878.0))[0]
print ind
print read[ind,1]
In [ ]:
indices = np.where(np.logical_and(np.in1d(read[:,1], tadds), np.in1d(read[:,2], tags)))[0]
indices.argsort
print sorted(indices)
print indices.shape
In [ ]:
print newtags
d = np.abs(read[:,0]-newtags[:,np.newaxis])
indices = np.where(np.any(np.isclose(d, 0.0), axis=0))[0]
#indices = np.where(np.isclose(read[:,0], 2.0))[0]
print indices
print indices.shape
print read[indices,0]
In [ ]:
aa = np.array([2,4,1,3])
bb = np.array([1,3,4,2])
for a in aa:
print np.where(bb==a)[0][0]
In [ ]:
#newtag = read[-1,0]
if newtag in read[:,0]:
indices = np.where(np.in1d(read[:,0], newtag))[0][0]
print indices
else:
print 'nothing'
In [ ]:
read[indices,:]
In [ ]:
den1 = read[indices,3]
tau1 = read[indices,4]
In [ ]:
tau1
In [ ]:
l
In [ ]: