Kelvin-Helmholtz Instability

In [1]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
In [2]:
nx=200
it=80
f1 =open('data/rho%05d.d'  %(it),'rb') 
f2 =open('data/vx%05d.d'   %(it),'rb')
f3 =open('data/vy%05d.d'   %(it),'rb')
f4 =open('data/vz%05d.d'   %(it),'rb') 
f5 =open('data/pr%05d.d'   %(it),'rb')
f6 =open('data/bx%05d.d'   %(it),'rb')
f7 =open('data/by%05d.d'   %(it),'rb') 
f8 =open('data/bz%05d.d'   %(it),'rb')    
f9 =open('data/beta%05d.d' %(it),'rb')
f10=open('data/eta%05d.d'  %(it),'rb')
f11=open('data/chi%05d.d'  %(it),'rb')
f12=open('data/jx%05d.d'   %(it),'rb')
f13=open('data/jy%05d.d'   %(it),'rb')
f14=open('data/jz%05d.d'   %(it),'rb')
f15=open('data/divu%05d.d' %(it),'rb')
f16=open('data/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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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='auto',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()

Sream lines

In [3]:
%%time
Lx=200.
Ly=500.
nx=200
ny=500
ni=1
nj=30
npops=2*ni*nj
nsteps=1500
len=.1
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]=100.0; fly[i*nj+j      ,0]=(Ly/2.0)/nj+Ly/nj*j
        flx[i*nj+j+ni*nj,0]=100.0; fly[i*nj+j+ni*nj,0]=(Ly/2.0)/nj+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*vx[iy+1,ix+1]+dx*(1-dy)*vx[iy,ix+1]+(1-dx)*dy*vx[iy+1,ix]+(1-dx)*(1-dy)*vx[iy,ix]
    sby=dx*dy*vy[iy+1,ix+1]+dx*(1-dy)*vy[iy,ix+1]+(1-dx)*dy*vy[iy+1,ix]+(1-dx)*(1-dy)*vy[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 here
    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

CPU times: user 249 ms, sys: 2 ms, total: 251 ms
Wall time: 251 ms
In [4]:
fig=plt.figure(figsize=(8,8))
#fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
#plt.imshow(rho,aspect='equal',origin='lower',cmap='jet');plt.colorbar(shrink=0.4)#;plt.title(r'$B$')
#plt.imshow(np.sqrt(by**2+bx**2),aspect='auto',origin='lower',cmap='bwr');plt.colorbar()#;plt.title(r'$B$')
plt.imshow(np.arctan(vy/vx)*180.0/np.pi,aspect='equal',origin='lower',cmap='bwr');plt.colorbar(shrink=0.8)#;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(50,150)
    #plt.ylim(0,100)

plt.xlabel(r'$x/\lambda_D$')
plt.ylabel(r'$y/\lambda_D$')
plt.show()

Animation

In [ ]:
!mkdir figures #if you have the directory, comment out this line.
In [ ]:
nx=200
for it in range(256,300):
    f  =open('data/rho%05d.d' %(it), 'rb')
    rho=np.fromfile(f, dtype='float64').reshape(-1,nx)
    
    fig=plt.figure(figsize=(8,8))
    #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.4)#;plt.title(r'$B$')
    plt.imshow(rho,aspect='auto',origin='lower',cmap='jet');plt.colorbar()#;plt.title(r'$B$')

    plt.savefig('figures/rho%03d.png' %(it), bbox_inches='tight',dpi=160)
    plt.close()
In [ ]:
!convert -delay 10 -loop 0 figures/rho*.png rho.gif

Animation is here!

In [ ]: