In [ ]:
import sys
import yt
from particles.particle_filters import *
import matplotlib.pyplot as plt
In [ ]:
#ts = yt.load('/d/d9/ychen/FLASH4/20180720_initM2/debug/GalaxyGroup_hdf5_chk_*')
#ts = yt.load('/d/d9/ychen/FLASH4/20180723_initBlockM7/data/GalaxyGroup_hdf5_plt_cnt_*')
#ts = yt.load('/d/d9/ychen/FLASH4/20180726_density_fix/debug/GalaxyGroup_hdf5_chk_0014.before.repeating.negative.dens')
#ds = yt.load('/d/d9/ychen/2018_production_runs/20180802_L438_rc10_beta07/debug/Group_L438_hdf5_plt_cnt_0924')
ds = yt.load('/d/d9/ychen/2018_production_runs/20180824_L438_rc30_beta07/data/Group_L438_hdf5_plt_cnt_0800')
print(ds)
In [ ]:
ds.print_stats()
In [ ]:
ds.index.level_stats['numgrids'][:12]
In [ ]:
print(ds)
#v, loc = ds.find_max('velocity_x')
#print('max velocity x:', v.in_units('c'), loc)
#v, loc = ds.find_max('velocity_y')
#print('max velocity y:', v.in_units('c'), loc)
v, loc = ds.find_max('velocity_z')
print('max velocity z:', v.in_units('c'), loc)\
v, loc = ds.find_max('alfven_speed')
print('max alfven speed:', v.in_units('c'), loc)
v, loc = ds.find_max('sound_speed')
print('max sound speed:', v.in_units('c'), loc)
v, loc = ds.find_max('velocity_magnitude')
print('max velocity:', v.in_units('c'), loc)
In [ ]:
#loc = ds.arr([1.688E+21, -1.194E+22, 3.694E+22], 'code_length')
#loc = ds.arr([1.194E+22, -1.812E+21, 2.681E+23], 'code_length')
loc = ds.arr([-1.438E+21, 3.938E+21, -4.062E+21], 'code_length')
#loc = ds.arr([-1.181E+22, 1.719E+22, -3.628E+23], 'code_length')
#loc = ds.arr([-1.938E+21, -2.312E+21, -3.994E+22], 'code_length')
v = ds.find_field_values_at_point('velocity_magnitude', loc).in_units('c')[0]
print('velocity_magnitude:', v.in_units('c'))
v = ds.find_field_values_at_point('alfven_speed', loc).in_units('c')[0]
print('alfven speed:', v.in_units('c'))
v = ds.find_field_values_at_point('sound_speed', loc).in_units('c')[0]
print('sound speed:', v.in_units('c'))
In [ ]:
fields = ['temperature', 'density', 'pressure','shks','velocity_z', 'velocity_y', 'jet ']
for axis in 'xyz':
plot = yt.SlicePlot(ds, axis, fields, center=loc, width=(10, 'kpc'), origin='native')
plot.annotate_marker(loc, plot_args={'color': 'k'})
#plot.set_zlim('jet ', 1E-10, 1E0)
plot.set_zlim('velocity_z', -6E9, 6E9)
plot.set_cmap('velocity_z', 'seismic')
plot.set_zlim('velocity_y', -6E8, 6E8)
plot.set_cmap('velocity_y', 'seismic')
#plot.set_cmap('magnetic_field_y', 'seismic')
#plot.set_zlim('magnetic_field_y', -2E-4, 2E-4)
#plot.set_zlim('alfven_speed', 1E7, 1E11)
#plot.set_zlim('shks', 1, 7)
plot.annotate_grids()#draw_ids=True)
#plot.annotate_particles()
plot.show()
In [ ]:
fields = ['temperature','density', 'pressure','shks','velocity_z', 'jet ']
#fields = ['eint']
for axis in 'xyz':
plot = yt.SlicePlot(ds, axis, fields, center=loc, width=(2, 'kpc'), origin='native')
plot.annotate_marker(loc, plot_args={'color': 'k'})
#plot.set_zlim('jet ', 1E-10, 1E0)
plot.set_zlim('velocity_z', -6E9, 6E9)
plot.set_cmap('velocity_z', 'seismic')
#plot.set_cmap('magnetic_field_y', 'seismic')
#plot.set_zlim('magnetic_field_y', -2E-4, 2E-4)
#plot.set_zlim('alfven_speed', 1E7, 1E11)
#plot.set_zlim('shks', 1, 7)
plot.annotate_grids()#draw_ids=True)
#plot.annotate_particles()
plot.show()
In [ ]:
In [ ]:
fields = ['velocity_z', 'eint']
#fields = ['eint']
for axis in 'xy':
plot = yt.SlicePlot(ds, axis, fields, center=loc, width=(1, 'kpc'), origin='native')
plot.annotate_marker(loc, plot_args={'color': 'k'})
#plot.set_zlim('jet ', 1E-10, 1E0)
plot.set_zlim('velocity_z', -6E9, 6E9)
plot.set_cmap('velocity_z', 'seismic')
plot.annotate_velocity()
plot.annotate_grids(draw_ids=True)
plot.show()
In [ ]:
ds2 = yt.load('/d/d9/ychen/2018_production_runs/20180801_L430_rc10_beta07/data/Group_L430_hdf5_plt_cnt_0410')
fields = ['velocity_z']
#fields = ['eint']
loc2 = loc+ds.arr([0, 1.225E+21, 0], 'code_length')
for axis in 'xy':
plot = yt.SlicePlot(ds2, axis, fields, center=loc2, width=(3, 'kpc'), origin='native')
plot.annotate_marker(loc2, plot_args={'color': 'k'})
#plot.set_zlim('jet ', 1E-10, 1E0)
plot.set_log('velocity_z', True, linthresh=1E8)
plot.set_zlim('velocity_z', -6E9, 6E9)
plot.set_cmap('velocity_z', 'seismic')
plot.annotate_velocity()
plot.annotate_grids()
plot.show()
In [ ]:
n = 1
p = ds.point(loc)
dx = p.fwidth[0,0]
print(dx.in_units('pc'))
le = loc - (n+0.5)*dx
print(loc, le)
In [ ]:
cg = ds.covering_grid(11, le, [n*2+1]*3)
print(cg.LeftEdge)
print(cg.center)
print('dvx = ', (cg['velocity_x'][2,1,1] - cg['velocity_x'][0,1,1]).in_units('c'))
print('dvy = ', (cg['velocity_y'][1,2,1] - cg['velocity_y'][1,0,1]).in_units('c'))
print('dvz = ', (cg['velocity_z'][1,1,2] - cg['velocity_z'][1,1,0]).in_units('c'))
print(cg['sound_speed'][1,1,1].in_units('c'))
In [ ]:
ds.max_level
In [ ]:
g = ds.index.grids[3563]
print(g.id)
i,j,k = 0,3,7
print((g['vely'][i,j+1,k]-g['vely'][i,j-1,k]).in_units('c'))
print(g['sound_speed'][i,j,k].in_units('c'))
In [ ]:
from yt.geometry.selection_routines import point_selector
from yt.data_objects.selection_data_containers import YTPoint
p = YTPoint(loc, ds)
In [ ]:
kpc = yt.units.kpc.in_units('cm')
fields = [
# 'thermal_energy',
# 'total_energy',
'eint',
# 'magnetic_energy',
# 'magnetic_field_magnitude',
'temperature',
'pressure',
'density',
'jet ',
'velocity_y'
# 'velocity_magnitude'
]
npoints = 4
lines = []
p = ds.point(loc)
dx = p.fwidth[0,0]
for axis in [0,1,2]:
posA = ds.arr(loc.copy(), 'code_length')
posA[axis] = posA[axis] - npoints*dx
posB = ds.arr(loc.copy(), 'code_length')
posB[axis] = posB[axis] + npoints*dx
lines.append(yt.LineBuffer(ds, posA, posB, 2*npoints+1, label='axis %i' % axis))
plot = yt.LinePlot.from_lines(ds, fields, lines)
plot.set_x_unit('kpc')
plot.annotate_legend(fields)
plot.show()
In [ ]:
np.log10(lines[0]['eint'])
In [ ]:
ds = ts
kpc = yt.units.kpc.in_units('cm')
fields = [
# 'thermal_energy',
'total_energy',
# 'magnetic_energy',
# 'magnetic_field_magnitude',
'temperature',
'pressure',
'density',
'jet ',
'velocity_y'
# 'velocity_magnitude'
]
npoints = 4
lines = []
p = ds.point(loc)
dx = p.fwidth[0,0]
for axis in [0,1,2]:
posA = ds.arr(loc.copy(), 'code_length')
posA[axis] = posA[axis] - npoints*dx
posB = ds.arr(loc.copy(), 'code_length')
posB[axis] = posB[axis] + npoints*dx
lines.append(yt.LineBuffer(ds, posA, posB, 2*npoints+1, label='axis %i' % axis))
plot = yt.LinePlot.from_lines(ds, fields, lines)
plot.set_x_unit('kpc')
plot.annotate_legend(fields)
plot.show()
In [ ]:
ad = ds.all_data()
plt.scatter(np.log10(ad['eint'][::100]), np.log10(np.abs(ad['z'][::100])), s=1, lw=0)
In [ ]:
grids = ds.index.grids
eint_maxes = []
eint_mins = []
for g in grids:
sys.stdout.write('\rLoading grid %6i /%6i)' % (g.id, len(grids)))
eint_maxes.append(np.max(g['eint'])/np.mean(g['eint']))
eint_mins.append(np.min(g['eint'])/np.mean(g['eint']))
In [ ]:
grids = ds.index.grids
std_maxes = []
std_mins = []
for g in grids:
sys.stdout.write('\rLoading grid %6i /%6i)' % (g.id, len(grids)))
std = np.std(np.log10(g['eint']))
mean = np.mean(np.log10(g['eint']))
std_maxes.append(np.max((np.log10(g['eint'])-mean)/std))
std_mins.append(np.min((np.log10(g['eint'])-mean)/std))
In [ ]:
mean
In [ ]:
std_maxes = np.array(std_maxes)
std_mins = np.array(std_mins)
print(sorted(std_maxes, reverse=True)[:20])
print(sorted(std_mins)[:20])
In [ ]:
plt.hist(std_mins)
plt.hist(std_maxes)
In [ ]:
eint_maxes = np.array(eint_maxes)
eint_mins = np.array(eint_mins)
In [ ]:
plt.hist(np.log10(eint_maxes))
plt.hist(np.log10(eint_mins))
In [ ]:
print(sorted(eint_maxes, reverse=True)[:10])
print(sorted(eint_mins)[:10])
In [ ]:
np.argwhere(np.array(eint_maxes) > 100)
In [ ]:
np.argwhere(eint_mins < 1/500)
In [ ]:
n = 45661
g = grids[n]
eint = g['eint']
print(eint)
print(g.id)
print(g.Level)
print('%e' % np.min(eint))
print('%e' % np.max(eint))
print('%e' % np.mean(eint))
print('%e' % np.std(eint))
print('%f' % (np.max(eint)/np.mean(eint)))
print('%f' % (np.min(eint)/np.mean(eint)))
print('%f' % ((np.max(eint)-np.mean(eint))/np.std(eint)))
print('%f' % ((np.min(eint)-np.mean(eint))/np.std(eint)))
In [ ]:
ind = (3,7,4)
gc = 2
g.retrieve_ghost_zones(gc, 'eint')['eint'][ind[0]:ind[0]+2*gc, ind[1]:ind[1]+2*gc, ind[2]:ind[2]+2*gc]
In [ ]:
for i in np.argsort(eint_mins)[:4]:
g = grids[i]
argmin = np.unravel_index(np.argmin(g['eint']), g.shape)
locmin = g.fcoords[np.argmin(g['eint'])]
print(i, '%e' % np.min(g['eint']), argmin, locmin)
plot = yt.SlicePlot(ds, 'x', 'eint', center=locmin, width=(1,'kpc'))
plot.annotate_grids()
plot.show()
In [ ]:
print(g['y'][7,7,7])
g.fcoords
In [ ]:
for i in np.argsort(eint_maxes)[-4:]:
g = grids[i]
argmax = np.argmax(g['eint'])
locmax = g.fcoords[argmax]
print(i, '%e' % np.max(g['eint']), np.unravel_index(argmax, g.shape), locmax)
plot = yt.SlicePlot(ds, 'x', 'eint', center=locmax, width=(2,'kpc'), origin='native')
plot.annotate_marker(locmax, plot_args={'color': 'k'})
plot.annotate_grids()
plot.show()
In [ ]:
plt.hist(np.log10(ad['eint']), bins=40)
In [ ]:
grid_levels = np.array([g.Level for g in grids])
In [ ]:
grid_levels[np.argwhere(eint_ratios > 1000)]
In [ ]:
plt.scatter(np.log10(ad['eint'][::100]), np.abs(ad['grid_level'][::100]), s=1, lw=0)
In [ ]:
plot = yt.SlicePlot(ds, 'x', 'pressure', width=(5, 'kpc'))
plot.set_zlim('pressure', 1E-11, 1E-9)
plot.annotate_grids()
plot.show()
In [ ]:
ds = ts[0]
plot = yt.SlicePlot(ds, 'x', 'pressure', width=(5,'kpc'))
plot.set_zlim('pressure', 1E-11, 1E-9)
plot.annotate_grids()
plot.show()
In [ ]:
for axis in 'xyz':
plot = yt.SlicePlot(ds, axis, 'shks', center=loc, width=(8, 'kpc'))
plot.annotate_marker(loc)
#plot.set_zlim('jet ', 1E-42, 1E0)
#plot.set_zlim('velocity_z', -6E9, 6E9)
#plot.set_cmap('velocity_z', 'seismic')
plot.annotate_grids()
plot.show()
In [ ]: