Blast shock wave + density flucutuations¶
In [2]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
In [3]:
nx=1024
it=210
f1 =open('data3/rho%05d.d' %(it),'rb')
f2 =open('data3/vx%05d.d' %(it),'rb')
f3 =open('data3/vy%05d.d' %(it),'rb')
f4 =open('data3/vz%05d.d' %(it),'rb')
f5 =open('data3/pr%05d.d' %(it),'rb')
f6 =open('data3/bx%05d.d' %(it),'rb')
f7 =open('data3/by%05d.d' %(it),'rb')
f8 =open('data3/bz%05d.d' %(it),'rb')
f9 =open('data3/beta%05d.d' %(it),'rb')
f10=open('data3/eta%05d.d' %(it),'rb')
f11=open('data3/chi%05d.d' %(it),'rb')
f12=open('data3/jx%05d.d' %(it),'rb')
f13=open('data3/jy%05d.d' %(it),'rb')
f14=open('data3/jz%05d.d' %(it),'rb')
f15=open('data3/divu%05d.d' %(it),'rb')
f16=open('data3/divb%05d.d' %(it),'rb')
rho =np.fromfile(f1, dtype='float64').reshape(-1,nx)
vx =np.fromfile(f2, dtype='float64').reshape(-1,nx)
vy =np.fromfile(f3, dtype='float64').reshape(-1,nx)
vz =np.fromfile(f4, dtype='float64').reshape(-1,nx)
pr =np.fromfile(f5, dtype='float64').reshape(-1,nx)
bx =np.fromfile(f6, dtype='float64').reshape(-1,nx)
by =np.fromfile(f7, dtype='float64').reshape(-1,nx)
bz =np.fromfile(f8, dtype='float64').reshape(-1,nx)
beta=np.fromfile(f9, dtype='float64').reshape(-1,nx)
eta =np.fromfile(f10,dtype='float64').reshape(-1,nx)
chi =np.fromfile(f11,dtype='float64').reshape(-1,nx)
jx =np.fromfile(f12,dtype='float64').reshape(-1,nx)
jy =np.fromfile(f13,dtype='float64').reshape(-1,nx)
jz =np.fromfile(f14,dtype='float64').reshape(-1,nx)
divu=np.fromfile(f15,dtype='float64').reshape(-1,nx)
divb=np.fromfile(f16,dtype='float64').reshape(-1,nx)
display.clear_output(True)
plt.figure(figsize=(20,12))
plt.subplots_adjust(wspace=0.2, hspace=0.2)
plt.subplot(3,5,1) ;plt.imshow(rho, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$\rho$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,2) ;plt.imshow(vx, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$v_x$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,3) ;plt.imshow(vy, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$v_y$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,4) ;plt.imshow(vz, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$v_z$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,5) ;plt.imshow(pr, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$Pr$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,6) ;plt.imshow(bx, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$B_x$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,7) ;plt.imshow(by, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$B_y$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,8) ;plt.imshow(bz, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$B_z$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,9) ;plt.imshow(beta,aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$\beta$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,10);plt.imshow(eta, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$\eta$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,11);plt.imshow(chi, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$\chi$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,12);plt.imshow(divb,aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$\nabla\cdot\bf{B}$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,13);plt.imshow(jx, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$j_x$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,14);plt.imshow(jy, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$j_y$'); plt.tick_params(labelbottom='off',labelleft='off')
plt.subplot(3,5,15);plt.imshow(jz, aspect='equal',origin='lower',cmap='jet');plt.colorbar();plt.title(r'$j_z$'); plt.tick_params(labelbottom='off',labelleft='off')
#plt.suptitle('it=%06d' %it)
plt.show()
In [4]:
fig=plt.figure(figsize=(12,12))
plt.imshow(np.sqrt(bx**2+by**2),cmap='hot',aspect='equal',origin='lower')
#plt.imshow(eta,cmap='Blues',aspect='equal')
#plt.contour(jz)
#plt.imshow(jz,cmap='bwr',aspect='equal',clim=(-80,80))
#plt.imshow(np.arctan(by/bx)*180/np.pi,aspect='equal',origin='lower',cmap='bwr')
#plt.imshow(rho,aspect='equal',cmap='Blues')
#plt.clim(0.,60)
plt.colorbar(shrink=0.8)
plt.show()
Time evolution of the magnetic field¶
In [5]:
b2hst=np.zeros((200))
nx=1024
for it in range(200):
f2 =open('data3/bx%05d.d' %(it),'rb')
f3 =open('data3/by%05d.d' %(it),'rb')
bx =np.fromfile(f2, dtype='float64').reshape(-1,nx)
by =np.fromfile(f3, dtype='float64').reshape(-1,nx)
b2hst[it]=np.mean(np.sqrt(by**2+bx**2))
plt.plot(b2hst)
plt.show()
Sream lines of the magnetic field¶
In [6]:
%%time
nx=1024
it=100
f6 =open('data3/bx%05d.d' %(it),'rb')
f7 =open('data3/by%05d.d' %(it),'rb')
bx =np.fromfile(f6, dtype='float64').reshape(-1,nx)
by =np.fromfile(f7, dtype='float64').reshape(-1,nx)
Lx=1024
Ly=1024
nx=1024
ny=1024
ni=1
nj=80
npops=2*ni*nj
nsteps=40000
len=.05
flx=np.zeros((npops,nsteps))
fly=np.zeros((npops,nsteps))
for i in range(ni):
for j in range(nj):
flx[i*nj+j ,0]=56; fly[i*nj+j ,0]=Ly/nj*(j)
flx[i*nj+j+ni*nj,0]=56; fly[i*nj+j+ni*nj,0]=Ly/nj*(j)
sgnx=np.r_[np.ones(npops//2), -np.ones(npops//2)]
sgny=np.r_[np.ones(npops//2), -np.ones(npops//2)]
for i in range(nsteps-1):
dx=flx[:,i]-np.int_(flx[:,i])
dy=fly[:,i]-np.int_(fly[:,i])
ix=np.int_(flx[:,i])
iy=np.int_(fly[:,i])
sbx=dx*dy*bx[iy+1,ix+1]+dx*(1-dy)*bx[iy,ix+1]+(1-dx)*dy*bx[iy+1,ix]+(1-dx)*(1-dy)*bx[iy,ix]
sby=dx*dy*by[iy+1,ix+1]+dx*(1-dy)*by[iy,ix+1]+(1-dx)*dy*by[iy+1,ix]+(1-dx)*(1-dy)*by[iy,ix]#dy*by[iy+1,ix+1]+(1-dy)*by[iy,ix]
newbx=sbx/np.sqrt(sbx**2+sby**2) #normalized
newby=sby/np.sqrt(sbx**2+sby**2) #normalized
flx[:,i+1]=flx[:,i]+sgnx*len*newbx
fly[:,i+1]=fly[:,i]+sgny*len*newby
#if it hits the boundary, it stays there
sgnx[fly[:,i+1]>=Ly-2]=0.0;sgnx[flx[:,i+1]>=Lx-2]=0.0
sgny[fly[:,i+1]>=Ly-2]=0.0;sgny[flx[:,i+1]>=Lx-2]=0.0
sgnx[fly[:,i+1]<=2]=0.0; sgnx[flx[:,i+1]<=2]=0.0
sgny[fly[:,i+1]<=2]=0.0; sgny[flx[:,i+1]<=2]=0.0
In [9]:
fig=plt.figure(figsize=(12,12))
#fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
#plt.imshow(np.arctan(by/bx)*180/np.pi,aspect='equal',origin='lower',cmap='bwr');plt.colorbar(shrink=0.8)#;plt.title(r'$B$')
#plt.imshow(np.sqrt(by**2+bx**2),aspect='equal',origin='lower',cmap='hot');plt.colorbar()#;plt.title(r'$B$')
for i in range(npops):
plt.plot(flx[i,:],fly[i,:],'k-',lw=1)
plt.plot(flx[i,0],fly[i,0],'r.',lw=1)
#plt.xlim(100,400)
#plt.ylim(100,400)
plt.xlabel(r'$x/\lambda_D$')
plt.ylabel(r'$y/\lambda_D$')
plt.show()
Animation¶
In [3]:
!mkdir figures3 #if you have the directory, comment out this line.
In [11]:
nx=1024
ny=1024
for it in range(1):
f1 =open('data3/rho%05d.d' %(it),'rb')
f2 =open('data3/bx%05d.d' %(it),'rb')
f3 =open('data3/by%05d.d' %(it),'rb')
f4 =open('data3/vx%05d.d' %(it),'rb')
f5 =open('data3/vy%05d.d' %(it),'rb')
f6 =open('data3/pr%05d.d' %(it),'rb')
rho =np.fromfile(f1, dtype='float64').reshape(-1,nx)
bx =np.fromfile(f2, dtype='float64').reshape(-1,nx)
by =np.fromfile(f3, dtype='float64').reshape(-1,nx)
vx =np.fromfile(f4, dtype='float64').reshape(-1,nx)
vy =np.fromfile(f5, dtype='float64').reshape(-1,nx)
pr =np.fromfile(f6, dtype='float64').reshape(-1,nx)
fig,((ax1,ax2),(ax3,ax4)) = plt.subplots(2, 2,figsize=(8,8))
fig.subplots_adjust(hspace=0,wspace=-0.05)
b2=np.sqrt(by**2+bx**2)
v2=np.sqrt(vy**2+vx**2)
ax1.imshow(b2[nx//2:nx,0:nx//2] ,cmap='hot' ,origin='lower')
ax2.imshow(pr[nx//2:nx,nx//2:nx],cmap='Blues' ,origin='lower')
ax3.imshow(v2[0:nx//2 ,0:nx//2] ,cmap='terrain_r' ,origin='lower')
ax4.imshow(rho[0:nx//2,nx//2:nx],cmap='gist_stern',origin='lower')
ax1.set_aspect('equal');ax2.set_aspect('equal');ax3.set_aspect('equal');ax4.set_aspect('equal')
ax1.axis('off'); ax2.axis('off'); ax3.axis('off'); ax4.axis('off')
ax1.text(0,ny//2-64 ,r'$|B|$' ,color='w',fontsize=14)
ax2.text(nx//2-64,ny//2-64,r'$P$' ,color='k',fontsize=14)
ax3.text(0,ny//2-64 ,r'$|v|$' ,color='k',fontsize=14)
ax4.text(ny//2-64,ny//2-64,r'$\rho$',color='w',fontsize=14)
plt.savefig('figures3/blastsk+turb%03d.png' %(it), bbox_inches='tight',dpi=160)
plt.close()
In [22]:
!convert -delay 10 -loop 0 figures3/blastsk+turb*.png blastsk+turb.gif
