In [ ]:
%matplotlib inline
import os
import sys
sys.path.append('/home/ychen/lib/util')
import util
import numpy as np
import logging
from matplotlib import pyplot as plt
import yt
#from yt.imods import *
logging.getLogger('yt').setLevel(logging.ERROR)
from tools import read_par
# Scan for files
dir = '/d/d9/ychen/FLASH4/MHD_Jet_3D/0429_test'
pars = read_par(dir)
def rescan(printlist=False):
files = util.scan_files(dir, '*hdf5_plt_cnt_[0-9][0-9][0-9][0-9]', printlist=printlist)
return files
files = rescan(True)
pos0 = [0.0, 0.0, 0.0]
pos1 = [0.0, \
pars['nozzleRadius']+0.5*pars['rFeatherOut'], \
0.0]
pos2 = [0.0, \
0.5*(pars['rFeatherIn']+pars['nozzleRadius']), \
pars['zTorInj']+0.5*pars['zFeather']]
#pos3 = [0.0, 1.5E21, 4.5E21]
pos3 = [0.0, \
0.5*(pars['rFeatherIn']+pars['nozzleRadius']), \
pars['zTorInj']+1.5*pars['zFeather']]
vecX = pars['nozzleVecX']
vecY = pars['nozzleVecY']
vecZ = pars['nozzleVecZ']
norm = np.sqrt(vecX**2+vecY**2+vecZ**2)
R = np.array([ [vecZ, vecY], [-vecY, vecZ] ])/norm
nR = pars['nozzleRadius']
nL = pars['nozzleHalfL']
HydroNozzleY = np.array([-nR, nR, nR, -nR])
HydroNozzleZ = np.array([ nL, nL, -nL, -nL])
nozzleCoords = [ np.dot(R, np.array([y,z])) for (y,z) in zip(HydroNozzleY,HydroNozzleZ) ]
In [ ]:
def sliceplot(field, zoom_fac):#, zlim1, zlim2):
alldata = ds.all_data()
mi, ma = alldata.quantities.extrema(field)
print field, mi, ma
slice = yt.SlicePlot(ds, 'x', field, center=center)#, width=(1.2E20,1.2E20))
if field in ['magnetic_field_x', 'magnetic_field_y', 'magnetic_field_z',\
'velocity_x', 'velocity_y', 'velocity_z', 'magx', 'magy', 'magz']:
slice.set_cmap(field, 'RdBu_r')
slice.set_log(field, False)
if field in ['magnetic_field_x']:
slice.set_zlim(field, -3.0E-5, 3.0E-5)
if field in ['velocity_y']:
slice.set_zlim(field, -3.0E8, 3.0E8)
slice.annotate_grids(edgecolors='black')
slice.zoom(zoom_fac)
slice.annotate_timestamp(0.15, 0.95, normalized=True, format="{time:3.1e} {units}")
slice.annotate_marker(pos0, marker='o', plot_args={'color':'white', 'edgecolor':'k'})
slice.annotate_marker(pos1, marker='o', plot_args={'color':'white', 'edgecolor':'k'})
slice.annotate_marker(pos2, marker='o', plot_args={'color':'white', 'edgecolor':'k'})
slice.annotate_marker(pos3, marker='o', plot_args={'color':'white', 'edgecolor':'k'})
for pos in nozzleCoords:
slice.annotate_marker(pos, marker='o', plot_args={'color':'white', 'edgecolor':'black'})
slice.show()
#slice.save(file.pathname)
return slice
#center = (-4.69E20,4.03e+21, 1.45e+22)
center = (0, 0, 0)
for file in files[-1:]:
print file.fullpath
ds = yt.load(file.fullpath)
print ds.current_time
slice = sliceplot('magnetic_field_x', 2)
slice = sliceplot('velocity_z', 2)
slice = sliceplot('velocity_y', 2)
In [ ]:
posA = [pos3[0], pos3[1], 2*pos3[2]]
posB = [pos3[0], pos3[1], -2*pos3[2]]
ray = ds.ray( posA, posB )
zs = np.linspace(posA[2], posB[2], len(ray['magnetic_field_x']))
plt.plot(zs, ray['magnetic_field_x'], color='red', label='Bx')
plt.ylim(-0.000008, 0.000008)
plt.legend(loc=2)
plt.twinx().plot(zs, ray['velocity_z'], color='blue', label='Vz')
plt.plot(zs, ray['velocity_y']*10, color='green', label='Vy')
plt.axhline(y=0)
plt.axvline(x=pos2[2])
plt.axvline(x=pos3[2])
plt.legend()
print pos2
In [ ]:
ts = []
P0s, P1s, P2s, P3s, Ps = [], [], [], [], []
rho0s = []
By0s, Bz0s, Bz1s, Bx2s, Bz2s, Bz3s = [], [], [], [], [], []
Ppols, Ptor2s, Ptor3s = [], [], []
for file in files[2:]:
#print file
ds = yt.load(file.fullpath)
ts.append(ds.current_time)
rho0 = ds.find_field_values_at_point('density', pos0)
P0 = ds.find_field_values_at_point('pressure', pos0)
P1 = ds.find_field_values_at_point('pressure', pos1)
P2 = ds.find_field_values_at_point('pressure', pos2)
P3 = ds.find_field_values_at_point('pressure', pos3)
rho0s.append(rho0)
P0s.append(P0)
P1s.append(P1)
P2s.append(P2)
P3s.append(P3)
Ps.append((P1+P2+P3)/3.0)
Bz0 = ds.find_field_values_at_point('magnetic_field_z', pos0)
By0 = ds.find_field_values_at_point('magnetic_field_y', pos0)
Bz1 = ds.find_field_values_at_point('magnetic_field_z', pos1)
Bx2 = ds.find_field_values_at_point('magnetic_field_x', pos2)
By2 = ds.find_field_values_at_point('magnetic_field_y', pos2)
Bz2 = ds.find_field_values_at_point('magnetic_field_z', pos2)
Bx3 = ds.find_field_values_at_point('magnetic_field_x', pos3)
By3 = ds.find_field_values_at_point('magnetic_field_y', pos3)
Bz3 = ds.find_field_values_at_point('magnetic_field_z', pos3)
Bz0s.append(Bz0)
Bx2s.append(Bx2)
Bz1s.append(Bz1)
Bz2s.append(Bz2)
Bz3s.append(Bz3)
Ppols.append((Bz0**2+By0**2)/8.0/np.pi)
Ptor2s.append((Bx2**2+By2**2)/8.0/np.pi)
Ptor3s.append((Bx3**2+By3**2)/8.0/np.pi)
In [ ]:
plt.plot(ts, P0s, '*-', label='nozzle')
#plt.plot(ts, Ps, 'o-', label='P')
plt.plot(ts, P1s, '.:', label='P1')
plt.plot(ts, P2s, '.:', label='P2')
plt.plot(ts, P3s, '.:', label='P3')
#plt.axhline(y=1.86E-6)
plt.legend(loc=0)
plt.xlabel('t')
plt.ylabel('pressure')
#plt.ylim(0,1.0E-6)
#plt.semilogy()
In [ ]:
plt.plot(ts, Ppols, '*-', label='Ppol_nozzle')
plt.plot(ts, Ptor2s, '.:', label='Ptor2')
plt.plot(ts, Ptor3s, '.:', label='Ptor3')
#plt.plot(ts, Bzs, 'o-', label='Bz')
#plt.axhline(y=1.86E-6)
plt.legend(loc=0)
plt.xlabel('t')
plt.ylabel('magnetic_field')
#plt.semilogy()
In [ ]:
Ppols = np.array(Ppols)
Ptor2s = np.array(Ptor2s)
Ptor3s = np.array(Ptor3s)
P0s = np.array(P0s)
beta = pars['sim_betaJet']
#beta = 20
plt.plot(ts, P0s, 'o-', label='Pgas')
plt.plot(ts, Ppols*beta, '*-', label='Ppol')
plt.plot(ts, Ptor2s*beta, '^-', label='Ptor2')
plt.plot(ts, Ptor3s*beta, 'v-', label='Ptor3')
plt.legend()
plt.show()
In [ ]:
#plt.plot(ts, P0s/(Ppols+Ptors))
plt.plot(ts, Ptor2s/Ppols, 'o-', label='tor2')
plt.plot(ts, Ptor3s/Ppols, 'o-', label='tor3')
pars['sim_helicityJet'] = 1
plt.axhline(y=pars['sim_helicityJet'])
plt.ylim(0.0,1.5)
plt.legend()
In [ ]:
# pressure
PI = np.pi
g = 1.6666
v = 3.0E9
r = 7.5E20
L = 1.0E45
sim_rhoAmbient = 1.67E-25
A = 2.0*PI*r*r
t0 = (A*v)**1.25*(sim_rhoAmbient*g/(g-1)/L)**0.75*0.227082*2.0
print 't0 = %e' % t0
time = np.linspace(0, r/v*2000, 1000)
tcut = np.array([t if t > t0 else t0 for t in time])
pressure = tcut**(-0.8)*0.305454*sim_rhoAmbient**0.6*((g-1)/g*L)**0.4
#plt.plot(time, P1s, '.:', label='P1')
#plt.plot(time, P2s, '.:', label='P2')
#plt.plot(time, P3s, '.:', label='P3')
plt.plot(time, pressure, '-')#, label='theory')
#plt.plot(ts, P0s, '*-')
#plt.axvline(x=t0, ls=':')
plt.legend(loc=0)
#plt.semilogy()
#plt.ylim(0.0, P0*20.0)
In [ ]:
# density
density = 2.0/v/v*(L/(A*v) - g/(g-1.0)*pressure)
plt.plot(time, density, '-', label='Jet density')
#plt.plot(ts, rho0s, '*-')
#plt.axhline(y=sim_rhoAmbient, label='Ambient density')
#plt.semilogy()
print '%e' % (L/(r*r*PI*2*v))
print '%e' % (g/(g-1.0)*pressure[0])
In [ ]:
# Mach number
Ms = v/np.sqrt(g*pressure/density)
plt.plot(time, Ms, '-', label='Mach Number')
print Ms[0]
In [ ]:
ekin = 0.5*density*v*v
eth = g/(g-1.0)*pressure
plt.plot(time, ekin, label='kinetic e')
plt.plot(time, eth, label='thermal e')
plt.semilogy()
plt.legend()
In [ ]:
t0 = (3.0/4.0/PI*A*v)**0.5*(A*v/L*g/(g-1.0)*sim_rhoAmbient/9.0)**0.75*2
print 't0 = %e' % t0
time = np.linspace(0, r/v*200, 1000)
tcut = np.array([t if t > t0 else t0 for t in time])
pressure = sim_rhoAmbient*(3.0/4.0/PI*A*v)**(2.0/3.0)*1.0/9.0*tcut**(-4.0/3.0)
plt.plot(time, pressure, '-')
In [ ]:
density = 2.0/v/v*(L/(A*v) - g/(g-1.0)*pressure)
plt.plot(time, density, '-', label='Jet density')
In [ ]:
# Mach number
Ms = v/np.sqrt(g*pressure/density)
plt.plot(time, Ms, '-', label='Mach Number')
print Ms[0]
In [ ]:
ekin = 0.5*density*v*v
eth = g/(g-1.0)*pressure
plt.plot(time, ekin, label='kinetic e')
plt.plot(time, eth, label='thermal e')
plt.semilogy()
plt.legend()
In [ ]:
PI = np.pi
g = 1.3333
v = 3.0E9
r = 7.5E20
bf =1.875E20
L = 1.0E45
M = np.linspace(0.5, 20, 100)
density = 0.5*L/PI/v**3/( 0.5*r*r*(1+1/M**2/(g-1)) + r*bf*(0.3125+1/M**2/(g-1)) + bf*bf*(0.06056+0.29736/M**2/(g-1)) )
pressure = v*v*density/M**2/g
plt.plot(M, pressure, label='pressure')
plt.semilogy()
plt.axhline(2.5E-10, ls=':', color='b')
plt.xlabel('M')
plt.ylabel('pressure')
plt.legend(loc=3)
plt.twinx()
plt.plot(M, density, label='density', color='r')
plt.semilogy()
plt.axhline(8.6E-26, ls=':', color='r')
plt.ylabel('density')
plt.legend(loc=1)
In [ ]:
kineticE = 0.5*density*v*v
thermalE = g/(g-1.0)*pressure
plt.plot(M, thermalE, label='thermal energy $E_{th}$')
plt.semilogy()
plt.xlabel('M')
plt.plot(M, kineticE, label='kinetic energy $E_k$')
plt.semilogy()
#plt.axvline(1.0/np.sqrt(g-1), ls=':', color='k')
plt.legend(loc=2)
plt.twinx()
plt.plot(M, kineticE/thermalE, color='r', label='$E_k/E_{th}$')
plt.plot(M, M, color='r', ls=':')
plt.legend(loc=1)
In [ ]:
sim_initMach = 5
sim_mach = 10
sim_tOn = 1.0E10
sim_duration = 1.58E15
time = np.linspace(0, 1.6E13, 999)
t1 = sim_duration/100.0
M = sim_initMach + (sim_mach-sim_initMach)*cos(PI*( 0.5*(time-sim_tOn-t1)/t1))
density = L/PI/v**3/( 0.5*r*r*(1+1/M**2/(g-1)) + r*bf*(0.3125+1/M**2/(g-1)) + bf*bf*(0.06056+0.29736/M**2/(g-1)) )
pressure = v*v*density/M**2/g
plt.plot(time, M, label='Mach', c='b')
plt.twinx()
plt.plot(time, pressure, c='r', label='pressure')
plt.twinx()
plt.plot(time, density, c='g', label='density')
In [ ]:
sim_velJet = 3.0E9
time = np.linspace(1.56E15, 1.60E15, 999)
#theta = 0.5*(time-sim_tOn-sim_duration+t1)/t1
#plt.plot(time, theta)
velocity = sim_velJet*cos(PI*( 0.5*(time-sim_tOn-sim_duration)/t1))
plt.plot(time, velocity)
In [ ]: