Magnetic Reconnection in a low-beta plasma

Reference is here (Zenitani2011PoP).

In [1]:
#Packages are defined here.
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
from subprocess import Popen

Simulation Parameters

In [2]:
nx=3200
ny=600
dx=0.1
dy=0.1
Lx=nx*dx  #->320
Ly=ny*dy  #->60

Snapshots

In [28]:
%%time
it=350
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()

CPU times: user 13.4 s, sys: 1.22 s, total: 14.6 s
Wall time: 14.6 s
In [29]:
fig=plt.figure(figsize=(15,10))
plt.imshow(vx,origin='lower',aspect='equal',cmap='Blues')
plt.colorbar(shrink=0.22)
#plt.ylim(0,300)
#plt.xlim(800,2400)
#plt.clim(-0.1,0.1)
plt.show()

Sream lines of the magnetic field

In [34]:
%%time
ni=1
nj=50
npops=2*ni*nj
nsteps=15000
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]=3100+nx/ni*i; fly[i*nj+j      ,0]=(ny/2.)/nj+ny/nj*j
        flx[i*nj+j+ni*nj,0]=3100+nx/ni*i; fly[i*nj+j+ni*nj,0]=(ny/2.)/nj+ny/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 here
    sgnx[fly[:,i+1]>=ny-2]=0.0;sgnx[flx[:,i+1]>=nx-2]=0.0 
    sgny[fly[:,i+1]>=ny-2]=0.0;sgny[flx[:,i+1]>=nx-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 2.78 s, sys: 1 ms, total: 2.78 s
Wall time: 2.78 s
In [38]:
fig=plt.figure(figsize=(8,8))
#fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
#plt.imshow(np.arctan2(by,bx)*180/np.pi,aspect='equal',origin='lower',cmap='bwr');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(vx,aspect='equal',origin='lower',cmap='Blues')
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(1600,3200)
    plt.ylim(0,300)

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 [24]:
%%bash
echo -e " 
import numpy as np
import matplotlib.pyplot as plt
nx=3200
for it in range(200,450):
    f2 =open('data/vx%05d.d'   %(it),'rb')
    #f3 =open('data/by%05d.d'   %(it),'rb')
    
    vx  =np.fromfile(f2, dtype='float64').reshape(-1,nx)
    #by  =np.fromfile(f3, dtype='float64').reshape(-1,nx)
    
    fig=plt.figure(figsize=(8,8))
    #fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
    plt.imshow(vx,aspect='equal',origin='lower',cmap='Blues');
    plt.colorbar(shrink=0.15)#;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.savefig('figures/vx%03d.png' %(it), bbox_inches='tight',dpi=160)
    plt.close()
" > animation.py
In [25]:
proc = Popen(["python", "animation.py"], shell=False)
In [26]:
!convert -delay 10 -loop 0 figures/vx*.png vx.gif

animation is here!

In [ ]: