In [1]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import datetime
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
print('Last modified =', "{0:%Y-%b-%d %H:%M:%S.%f}".format(datetime.datetime.now()))

Last modified = 2018-May-06 14:09:22.825334

pCANS 2D Magnetic reconnection

Document is here: http://www.astro.phys.s.chiba-u.ac.jp/pcans/em2d_sk.html
note: I modified fio.f90 for output fies (see the bottom)
note: particle data is output as momentum
Physical paramterers:

  • $\Delta x=\Delta y=\lambda_D=1$
  • $\Delta t=0.5\Delta x/c$
  • $m_i/m_e$=16
  • $\omega_{pe}/\Omega_e=2$
  • $T_e/T_i=2$

Simulation parameters

In [3]:
nx=320
ny=10
lx=320.0
ly=10.0
nt=30000
dt=0.5*0.05  #dt = 0.5*0.05(c/\omega_{pe})
intvl=500
nrank=64
nd=nt//intvl  #the number of data
rmass=16
dir='data' #the directory location of the data

2D plot of the electromagnetic fields

In [4]:
%%time
it=10000
nmv=['bx','by','bz','ex','ey','ez']
bx,by,bz,ex,ey,ez=np.zeros((6,ny*nrank,nx))
EBs=[bx,by,bz,ex,ey,ez]
for i,EB in enumerate(EBs):
    for irank in range(nrank):
        f=open('%s/%s%06d_rank=%03d.dat' %(dir,nmv[i],it,irank),'rb')
        d=np.fromfile(f,dtype='float64').reshape(-1,nx+3)
        EB[ny*irank:ny*(irank+1),:]=d[0:ny,0:nx]

CPU times: user 40 ms, sys: 49 ms, total: 89 ms
Wall time: 4.73 s
In [5]:
%%time
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(12,16))
plt.subplots_adjust(hspace=0.2,wspace=0.5)
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt))
                  
im1=ax1.imshow(bx,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im1,ax=ax1);ax1.set_title(r'$B_x$')
im2=ax2.imshow(by,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im2,ax=ax2);ax2.set_title(r'$B_y$')
im3=ax3.imshow(bz,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im3,ax=ax3);ax3.set_title(r'$B_z$')
im4=ax4.imshow(ex,aspect='auto',origin='lower',cmap='terrain_r',clim=(-.5,.5));plt.colorbar(im4,ax=ax4);ax4.set_title(r'$E_x$')
im5=ax5.imshow(ey,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im5,ax=ax5);ax5.set_title(r'$E_y$')
im6=ax6.imshow(ez,aspect='auto',origin='lower',cmap='terrain_r');plt.colorbar(im6,ax=ax6);ax6.set_title(r'$E_z$')
ax1.set_xlabel(r'$x/\lambda_D$');ax2.set_xlabel(r'$x/\lambda_D$');ax3.set_xlabel(r'$x/\lambda_D$');ax4.set_xlabel(r'$x/\lambda_D$');ax5.set_xlabel(r'$x/\lambda_D$');ax6.set_xlabel(r'$x/\lambda_D$')
ax1.set_ylabel(r'$y/\lambda_D$');ax2.set_ylabel(r'$y/\lambda_D$');ax3.set_ylabel(r'$y/\lambda_D$');ax4.set_ylabel(r'$y/\lambda_D$');ax5.set_ylabel(r'$y/\lambda_D$');ax6.set_ylabel(r'$y/\lambda_D$')

plt.show()

CPU times: user 5.17 s, sys: 174 ms, total: 5.34 s
Wall time: 5.35 s

Stream lines of magnetic fields

In [6]:
%%time
Lx=320
Ly=ly*nrank
ni=25
nj=2
npops=2*ni*nj
nsteps=1000
len=.5
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]=10+Lx/ni*(i); fly[i*nj+j      ,0]=160+320*j
        flx[i*nj+j+ni*nj,0]=10+Lx/ni*(i); fly[i*nj+j+ni*nj,0]=160+320*j
        #print(i,j,i*nj+j,i*nj+j+ni*nj)
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]>=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 130 ms, sys: 2 ms, total: 132 ms
Wall time: 130 ms
In [7]:
fig=plt.figure(figsize=(4,8))
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
plt.imshow(np.sqrt(bx**2+by**2),aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,ny*nrank]);plt.colorbar()#;plt.title(r'$B$')
for i in range(npops):
    plt.plot(flx[i,:],fly[i,:],'w-',lw=.5)
    plt.plot(flx[i,0],fly[i,0],'r.',lw=1)

    #plt.xlim(50,270)

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

2D snapshot of density, average velocity and pressure

In [8]:
%%time
it=10000
binx=160
biny=10
nmv=['up1','up2']
mass=[rmass,1]
flds=np.zeros((8,biny*nrank,binx))

for i in range(2):
    for irank in range(nrank):   

        f=open('%s/%s%06d_rank=%03d.dat' %(dir,nmv[i],it,irank),'rb')
        up =np.fromfile(f,dtype='float64').reshape(-1,5)
        gam=np.sqrt(1+up[:,2]**2+up[:,3]**2+up[:,4]**2)
        for j in range(2,5): up[:,j]=up[:,j]/gam

        den,xedges,yedges =np.histogram2d(up[:,0],up[:,1],bins=(binx,biny))
        avlx,xedges,yedges=np.histogram2d(up[:,0],up[:,1],bins=(binx,biny),weights=up[:,2]); avlx=avlx/den
        avly,xedges,yedges=np.histogram2d(up[:,0],up[:,1],bins=(binx,biny),weights=up[:,3]); avly=avly/den
        pr=0.5*((up[:,2]-avlx[np.int_((up[:,0]-2)/lx*binx),np.int_((up[:,1]-2-ly*irank)/(ny+2)*biny)])**2+(up[:,3]-avly[np.int_((up[:,0]-2)/lx*binx),np.int_((up[:,1]-2-ly*irank)/(ny+2)*biny)])**2+up[:,4]**2)/3.0*mass[i]
        pr ,xedges,yedges=np.histogram2d(up[:,0],up[:,1],bins=(binx,biny),weights=pr)      

        flds[i  ,biny*irank:biny*(irank+1),:]=den.T
        flds[i+2,biny*irank:biny*(irank+1),:]=avlx.T
        flds[i+4,biny*irank:biny*(irank+1),:]=avly.T
        flds[i+6,biny*irank:biny*(irank+1),:]=pr.T

Ni=flds[0];Ne=flds[1];Vix=flds[2];Vex=flds[3];Viy=flds[4];Vey=flds[5];Pi=flds[6];Pe=flds[7]

CPU times: user 4.86 s, sys: 131 ms, total: 4.99 s
Wall time: 5.04 s
In [9]:
%%time
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(12,16))
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16,y=0.92)
plt.subplots_adjust(wspace=0.6)

im1=ax1.imshow(Ni ,aspect='auto',origin='lower',cmap='jet',extent=[0,nx,0,ny*nrank]);plt.colorbar(im1,ax=ax1);ax1.set_title(r'$N_i$')
im2=ax2.imshow(Viy,aspect='auto',origin='lower',cmap='jet',extent=[0,nx,0,ny*nrank]);plt.colorbar(im2,ax=ax2);ax2.set_title(r'$V_{iy}$')
im3=ax3.imshow(Pi ,aspect='auto',origin='lower',cmap='jet',extent=[0,nx,0,ny*nrank]);plt.colorbar(im3,ax=ax3);ax3.set_title(r'$P_i$')
im4=ax4.imshow(Ne ,aspect='auto',origin='lower',cmap='jet',extent=[0,nx,0,ny*nrank]);plt.colorbar(im4,ax=ax4);ax4.set_title(r'$N_e$')
im5=ax5.imshow(Vey,aspect='auto',origin='lower',cmap='jet',extent=[0,nx,0,ny*nrank]);plt.colorbar(im5,ax=ax5);ax5.set_title(r'$V_{ey}$')
im6=ax6.imshow(Pe ,aspect='auto',origin='lower',cmap='jet',extent=[0,nx,0,ny]);plt.colorbar(im6,ax=ax6);ax6.set_title(r'$P_e$')
ax1.set_xlabel(r'$x/\lambda_D$');ax2.set_xlabel(r'$x/\lambda_D$');ax3.set_xlabel(r'$x/\lambda_D$');ax4.set_xlabel(r'$x/\lambda_D$');ax5.set_xlabel(r'$x/\lambda_D$');ax6.set_xlabel(r'$x/\lambda_D$')
ax1.set_ylabel(r'$y/\lambda_D$');ax2.set_ylabel(r'$y/\lambda_D$');ax3.set_ylabel(r'$y/\lambda_D$');ax4.set_ylabel(r'$y/\lambda_D$');ax5.set_ylabel(r'$y/\lambda_D$');ax6.set_ylabel(r'$y/\lambda_D$')

plt.show()

CPU times: user 5.26 s, sys: 28 ms, total: 5.29 s
Wall time: 5.29 s
In [10]:
#integrating over the x-axis
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(15,5))
plt.subplots_adjust(hspace=0.4,wspace=0.3)
im1=ax1.plot(np.linspace(0,ly*nrank,biny*nrank),np.mean(Ni ,axis=1));ax1.set_xlabel(r'$y/\lambda_D$')     ;ax1.set_ylabel(r'$N_i$')
im2=ax2.plot(np.linspace(0,ly*nrank,biny*nrank),np.mean(Viy,axis=1));ax2.set_xlabel(r'$y/\lambda_D$');ax2.set_ylabel(r'$V_{iy}/c$')
im3=ax3.plot(np.linspace(0,ly*nrank,biny*nrank),np.mean(Pi ,axis=1));ax3.set_xlabel(r'$y/\lambda_D$')     ;ax3.set_ylabel(r'$P_i$')
im4=ax4.plot(np.linspace(0,ly*nrank,biny*nrank),np.mean(Ne ,axis=1));ax4.set_xlabel(r'$y/\lambda_D$')     ;ax4.set_ylabel(r'$N_e$')
im5=ax5.plot(np.linspace(0,ly*nrank,biny*nrank),np.mean(Vey,axis=1));ax5.set_xlabel(r'$y/\lambda_D$');ax5.set_ylabel(r'$V_{ey}/c$')
im6=ax6.plot(np.linspace(0,ly*nrank,biny*nrank),np.mean(Pe ,axis=1));ax6.set_xlabel(r'$y/\lambda_D$')     ;ax6.set_ylabel(r'$P_e$')

plt.show()

Snapshot of the ideal MHD approximation

In [11]:
plt.figure(figsize=(4,8))
plt.imshow(ez[:,::2]+(Vix*by[:,::2]-Viy*bx[:,::2]) ,aspect='auto',origin='lower',cmap='bwr',extent=[0,nx,0,ny*nrank]);
plt.title(r'$[{\bf E}+{\bf V}_i\times{\bf B}]_z$')
plt.clim(-.5,.5)
plt.xlabel(r'$x/\lambda_D$')
plt.ylabel(r'$y/\lambda_D$');
plt.colorbar()
plt.show()

Animation

In [ ]:
!mkdir figures #if you have the directory, comment out this line.
In [ ]:
%%time
nmv=['ex','ey','ez','bx','by','bz']
EB=np.zeros((6,ny*nrank,nx))
for iit in range(nd):
    it=(iit+1)*intvl
    for i in range(6):
        for irank in range(nrank):
            f=open('%s/%s%06d_rank=%03d.dat' %(dir,nmv[i],it,irank),'rb')
            d=np.fromfile(f,dtype='float64').reshape(-1,nx+3)
            EB[i,ny*irank:ny*(irank+1),:]=d[0:ny,0:nx]
    
    Lx=320
    Ly=ly*nrank
    ni=20
    nj=2
    npops=2*ni*nj
    nsteps=1000
    len=.5
    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]=10+Lx/ni*(i); fly[i*nj+j      ,0]=160+320*j
            flx[i*nj+j+ni*nj,0]=10+Lx/ni*(i); fly[i*nj+j+ni*nj,0]=160+320*j
            #print(i,j,i*nj+j,i*nj+j+ni*nj)
    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*EB[3,iy+1,ix+1]+(1-dx)*EB[3,iy,ix]
        sby=dy*EB[4,iy+1,ix+1]+(1-dy)*EB[4,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

        #flx[fly[:,it]>=ly-3,it+1]=flx[fly[:,it]>=ly-3,it] #if it hits the boundary, it stays there.    
        #fly[fly[:,it]>=ly-3,it+1]=fly[fly[:,it]>=ly-3,it] #if it hits the boundary, it stays there.

        #if it hits the boundary , it turns opposite
        sgnx[fly[:,i+1]>=Ly-2]=-sgnx[fly[:,i+1]>=Ly-2] 
        sgny[fly[:,i+1]>=Ly-2]=-sgny[fly[:,i+1]>=Ly-2] 
        sgnx[fly[:,i+1]<=2]   =-sgnx[fly[:,i+1]<=2]    
        sgny[fly[:,i+1]<=2]   =-sgny[fly[:,i+1]<=2]    
        flx[fly[:,i+1]>=Ly-2,i+1]=flx[fly[:,i+1]>=Ly-2,i]
        fly[fly[:,i+1]>=Ly-2,i+1]=fly[fly[:,i+1]>=Ly-2,i]
        flx[fly[:,i+1]<=2   ,i+1]=flx[fly[:,i+1]<=2   ,i]
        fly[fly[:,i+1]<=2   ,i+1]=fly[fly[:,i+1]<=2   ,i]    

    fig=plt.figure(figsize=(8,16))
    plt.imshow(EB[5,:,:],aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,ny*nrank])
    plt.title(r"$t\omega_{pe}=%.2f$" %(it*dt))

    for i in range(npops):
        plt.plot(flx[i,:],fly[i,:],'k-',lw=.5)

    plt.xlabel(r'$x/\lambda_D$')
    plt.ylabel(r'$y/\lambda_D$')
    plt.colorbar(label=r'$|B|$')

    fig.savefig('figures/B+strlin%03d.png' %(iit), bbox_inches='tight')
    ##print(iit)
    plt.close()
In [ ]:
!convert -delay 15 -loop 0 figures/B+strlin*.png  B+strlin.gif

Animation is here!

fio.f90

Part of fio.f90 is modified as follows

write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'bx',it2,'_rank=',nrank,'.dat'
    open(110+nrank,file=fname,form='unformatted',access="stream")
    do i=nys-1,nye+1
      write(110+nrank) (uf(1,j,i), j =nxs-1,nxe+1) 
    end do
    close(110+nrank)

    write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'by',it2,'_rank=',nrank,'.dat'
    open(110+nrank,file=fname,form='unformatted',access="stream")
    do i=nys-1,nye+1
      write(110+nrank) (uf(2,j,i), j=nxs-1,nxe+1) 
    end do
    close(110+nrank)

    write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'bz',it2,'_rank=',nrank,'.dat'
    open(110+nrank,file=fname,form='unformatted',access="stream")
    do i=nys-1,nye+1
      write(110+nrank) (uf(3,j,i), j=nxs-1,nxe+1) 
    end do
    close(110+nrank)

    write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'ex',it2,'_rank=',nrank,'.dat'
    open(110+nrank,file=fname,form='unformatted',access="stream")
    do i=nys-1,nye+1
      write(110+nrank) (uf(4,j,i), j=nxs-1,nxe+1) 
    end do
    close(110+nrank)

    write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'ey',it2,'_rank=',nrank,'.dat'
    open(110+nrank,file=fname,form='unformatted',access="stream")
    do i=nys-1,nye+1
      write(110+nrank) (uf(5,j,i), j=nxs-1,nxe+1) 
    end do 
    close(110+nrank)

    write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'ez',it2,'_rank=',nrank,'.dat'
    open(110+nrank,file=fname,form='unformatted',access="stream")
    do i=nys-1,nye+1
      write(110+nrank) (uf(6,j,i), j=nxs-1,nxe+1) 
    end do
    close(110+nrank)

    if(mod(it2,5000)==0) then
      write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'up1',it2,'_rank=',nrank,'.dat'
      open(110+nrank,file=fname,form='unformatted',access="stream")
      do j=nys,nye
        do ii=1,np2(j,1)
         write(110+nrank) (up(i,ii,j,1), i=1,5) 
        end do
      end do
      close(110+nrank)

      write(fname,'(a,a,i6.6,a,i3.3,a)') trim(dir),'up2',it2,'_rank=',nrank,'.dat'
      open(110+nrank,file=fname,form='unformatted',access="stream")
      do j=nys,nye
        do ii=1,np2(j,2)
         write(110+nrank) (up(i,ii,j,2), i=1,5) 
        end do
      end do
      close(110+nrank)
    endif