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:06:30.606707

pCANS 2D collisionless shock wave

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$=25
  • $\omega_{pe}/\Omega_e=10$
  • ${\bf B}_0=B_0{\hat y}$
  • $M_A=14$
  • $\Theta_{Bn}=90^\circ$
  • $\beta_i=\beta_e=0.5$

Simulation parameters

In [3]:
nx=4096
ny=8
lx=4096
ly=8.0
nt=60000
dt=0.5*0.05  #dt = 0.5*0.05(c/\omega_{pe})
intvl=500
nrank=32
nd=nt//intvl  #the number of data
rmass=25
dir='data' #the directory location of the data

2D plot of the electromagnetic fields

In [4]:
%%time
it=60000
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+2)
        EB[ny*irank:ny*(irank+1),:]=d[0:ny,0:nx]

CPU times: user 43 ms, sys: 57 ms, total: 100 ms
Wall time: 13.2 s
In [5]:
%%time
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(16,8))
plt.subplots_adjust(hspace=0.5,wspace=0.3)
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
                  
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=(-1,1));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.25 s, sys: 534 ms, total: 5.78 s
Wall time: 5.78 s
In [6]:
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(15,5))
plt.subplots_adjust(hspace=0.4)
im1=ax1.plot(np.mean(bx,axis=0));ax1.set_title('Bx')
im2=ax2.plot(np.mean(by,axis=0));ax2.set_title('By')
im3=ax3.plot(np.mean(bz,axis=0));ax3.set_title('Bz')
im4=ax4.plot(np.mean(ex,axis=0));ax4.set_title('Ex')
im5=ax5.plot(np.mean(ey,axis=0));ax5.set_title('Ey')
im6=ax6.plot(np.mean(ez,axis=0));ax6.set_title('Ez')
ax4.set_ylim(-1,1)

plt.show()

Stream lines of magnetic fields

In [13]:
%%time
Lx=4096
Ly=ly*nrank
npops=100
nsteps=5000
len=.1
flx=np.zeros((npops,nsteps))
fly=np.zeros((npops,nsteps))
flx[:,0]=np.linspace(1000,4000,npops)#np.random.rand(npops)*lx
fly[:,0]=3                           #np.random.rand(npops)*ly
sgnx=np.ones(npops)
sgny=np.ones(npops)
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 603 ms, sys: 3 ms, total: 606 ms
Wall time: 605 ms
In [14]:
fig=plt.figure()
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,:],'k-',lw=.5)
    plt.plot(flx[i,0],fly[i,0],'r.',lw=1)

    plt.xlim(1000,3500)

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

Stacked plot of the magnetic field $\left|{\bf B}\right|$

In [9]:
%%time

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 it in range(nd):
    iit=intvl*(it+1)
    for i,EB in enumerate(EBs):
        for irank in range(nrank):
            f=open('%s/%s%06d_rank=%03d.dat' %(dir,nmv[i],iit,irank),'rb')
            d=np.fromfile(f,dtype='float64').reshape(-1,nx+2)
            EB[ny*irank:ny*(irank+1),:]=d[0:ny,0:nx]
            
    plt.plot(np.mean(np.sqrt(bx**2+by**2+bz**2),axis=0)+it,'k',lw=0.5)
    plt.xlim(0,4000)

plt.show()

CPU times: user 11.7 s, sys: 6.47 s, total: 18.2 s
Wall time: 26min 4s

Particle data

In [10]:
%%time
# load particle data
it=60000
for irank in range(1):
    f1=open('%s/up1%06d_rank=%03d.dat' %(dir,it,irank),'rb')
    f2=open('%s/up2%06d_rank=%03d.dat' %(dir,it,irank),'rb')
    
    up1=np.fromfile(f1,dtype='float64').reshape(-1,5)
    up2=np.fromfile(f2,dtype='float64').reshape(-1,5)
    gam1=np.sqrt(1+up1[:,2]**2+up1[:,3]**2+up1[:,4]**2)
    gam2=np.sqrt(1+up2[:,2]**2+up2[:,3]**2+up2[:,4]**2)
    for i in range(2,5): up1[:,i]=up1[:,i]/gam1  #transform to velocity
    for i in range(2,5): up2[:,i]=up2[:,i]/gam2  #transform to velocity

CPU times: user 521 ms, sys: 198 ms, total: 719 ms
Wall time: 739 ms

Phase-space plot

In [11]:
%%time
ph=np.zeros((8,500,200))
gam1=1/np.sqrt(1-up1[:,2]**2-up1[:,3]**2-up1[:,4]**2)-1
gam2=1/np.sqrt(1-up2[:,2]**2-up2[:,3]**2-up2[:,4]**2)-1
ph[0,:,:],xedges,yedges=np.histogram2d(up1[:,0],up1[:,2],bins=(500,200),range=[[0,nx],[-1,1]])
ph[1,:,:],xedges,yedges=np.histogram2d(up1[:,0],up1[:,3],bins=(500,200),range=[[0,nx],[-1,1]])
ph[2,:,:],xedges,yedges=np.histogram2d(up1[:,0],up1[:,4],bins=(500,200),range=[[0,nx],[-1,1]])
ph[3,:,:],xedges,yedges=np.histogram2d(up2[:,0],up2[:,2],bins=(500,200),range=[[0,nx],[-1,1]])
ph[4,:,:],xedges,yedges=np.histogram2d(up2[:,0],up2[:,3],bins=(500,200),range=[[0,nx],[-1,1]])
ph[5,:,:],xedges,yedges=np.histogram2d(up2[:,0],up2[:,4],bins=(500,200),range=[[0,nx],[-1,1]])
ph[6,:,:],xedges,yedges=np.histogram2d(up1[:,0],gam1    ,bins=(500,200),range=[[0,nx],[0.01,10]])
ph[7,:,:],xedges,yedges=np.histogram2d(up2[:,0],gam2    ,bins=(500,200),range=[[0,nx],[0.01,10]])

fig=plt.figure(figsize=(15,6))
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
ax1 = fig.add_subplot(421)
ax2 = fig.add_subplot(423, sharex=ax1, sharey=ax1)
ax3 = fig.add_subplot(425, sharex=ax1, sharey=ax1)
ax4 = fig.add_subplot(422, sharex=ax1, sharey=ax1)
ax5 = fig.add_subplot(424, sharex=ax1, sharey=ax1)
ax6 = fig.add_subplot(426, sharex=ax1, sharey=ax1)
ax7 = fig.add_subplot(427)
ax8 = fig.add_subplot(428)

im1=ax1.imshow(ph[0,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,-1,1])
im2=ax2.imshow(ph[1,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,-1,1])
im3=ax3.imshow(ph[2,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,-1,1])
im4=ax4.imshow(ph[3,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,-1,1])
im5=ax5.imshow(ph[4,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,-1,1])
im6=ax6.imshow(ph[5,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,-1,1])
im7=ax7.imshow(ph[6,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,0.01,10])
im8=ax8.imshow(ph[7,:,:].T+0.01,aspect='auto',origin='lower',cmap='terrain_r', norm=mpl.colors.LogNorm(),interpolation="gaussian",extent=[0,nx,0.01,10])

ax1.set_ylabel(r'$v_x/c$');ax2.set_ylabel(r'$v_y/c$');ax3.set_ylabel(r'$v_z/c$');ax4.set_ylabel(r'$v_x/c$');ax5.set_ylabel(r'$v_y/c$');ax6.set_ylabel(r'$v_z/c$')
ax7.set_ylabel(r'$\epsilon$');ax8.set_ylabel(r'$\epsilon$')
ax7.set_xlabel(r'$x/(c/\omega_{e})$');ax8.set_xlabel(r'$x/(c/\omega_{e})$')

ax1.set_xticklabels(())
ax7.set_yscale('log')
ax8.set_yscale('log')
ax1.locator_params(nbins=4)

plt.colorbar(im1,ax=ax1)
plt.colorbar(im2,ax=ax2)
plt.colorbar(im3,ax=ax3)
plt.colorbar(im4,ax=ax4)
plt.colorbar(im5,ax=ax5)
plt.colorbar(im6,ax=ax6)
plt.colorbar(im7,ax=ax7)
plt.colorbar(im8,ax=ax8)

#fig.savefig('1sk30eleionph.pdf', bbox_inches='tight',dpi=200)
plt.show()

CPU times: user 9.24 s, sys: 301 ms, total: 9.54 s
Wall time: 9.54 s
Parser   : 162 ms

Velocity distribution function $(v_x, N(v_x))$

In [12]:
bs=1500
xmin=1000
xmax=2000
plt.figure(figsize=(15,5))
plt.subplot(121)
plt.plot(up2[:,0],up2[:,2],',',color='C1')
plt.plot(up1[:,0],up1[:,2],',',color='C0')
plt.fill([xmin,xmax,xmax,xmin],[-1,-1,1,1],color="k",alpha=0.4)
plt.ylim(-1,1)
plt.subplot(122)
plt.hist(up1[(up1[:,0]>=xmin)&(up1[:,0]<xmax),2],bins=bs,range=(-0.25,0.25),histtype='step',color='C0')
plt.hist(up2[(up2[:,0]>=xmin)&(up2[:,0]<xmax),2],bins=bs,range=(-1   ,   1),histtype='step',color='C1')
plt.xlabel(r'$v_x/c$')
plt.ylabel(r'${\rm number\ of\ particles}$')
plt.show()

2D snapshot of density, average velocity and pressure

In [13]:
%%time
it=60000
binx=200
biny=10
nmv=['up1','up2']
mass=[rmass,1]
flds=np.zeros((6,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))
        avl,xedges,yedges=np.histogram2d(up[:,0],up[:,1],bins=(binx,biny),weights=up[:,2]); avl=avl/den
        pr=0.5*((up[:,2]-avl[np.int_((up[:,0]-2)/lx*binx),np.int_((up[:,1]-2-ly*irank)/ly/nrank*biny)])**2+up[:,3]**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),:]=avl.T
        flds[i+4,biny*irank:biny*(irank+1),:]=pr.T

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

CPU times: user 1min 40s, sys: 10.9 s, total: 1min 51s
Wall time: 1min 51s
In [14]:
%%time
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(16,8))
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)
plt.subplots_adjust(hspace=0.4,wspace=0.3)

im1=ax1.imshow(Ni ,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,(ny+2)*nrank]);plt.colorbar(im1,ax=ax1);ax1.set_title(r'$N_i$')
im2=ax2.imshow(Vix,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,(ny+2)*nrank]);plt.colorbar(im2,ax=ax2);ax2.set_title(r'$V_i$')
im3=ax3.imshow(Pi ,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,(ny+2)*nrank]);plt.colorbar(im3,ax=ax3);ax3.set_title(r'$P_i$')
im4=ax4.imshow(Ne ,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,(ny+2)*nrank]);plt.colorbar(im4,ax=ax4);ax4.set_title(r'$N_e$')
im5=ax5.imshow(Vex,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,nx,0,(ny+2)*nrank]);plt.colorbar(im5,ax=ax5);ax5.set_title(r'$V_e$')
im6=ax6.imshow(Pe ,aspect='auto',origin='lower',cmap='terrain_r',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 3.95 s, sys: 10 ms, total: 3.96 s
Wall time: 3.96 s
In [15]:
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,lx,binx),np.mean(Ni ,axis=0));ax1.set_ylabel(r'$N_i$')
im2=ax2.plot(np.linspace(0,lx,binx),np.mean(Vix,axis=0));ax2.set_ylabel(r'$V_{ix}/c$')
im3=ax3.plot(np.linspace(0,lx,binx),np.mean(Pi ,axis=0));ax3.set_ylabel(r'$P_i$')
im4=ax4.plot(np.linspace(0,lx,binx),np.mean(Ne ,axis=0));ax4.set_ylabel(r'$N_e$')
im5=ax5.plot(np.linspace(0,lx,binx),np.mean(Vex,axis=0));ax5.set_ylabel(r'$V_{ex}/c$')
im6=ax6.plot(np.linspace(0,lx,binx),np.mean(Pe ,axis=0));ax6.set_ylabel(r'$P_e$')

plt.show()

Energy spectrum

In [16]:
bins=200
xmin=2000
xmax=4000
enmin=1e-6
enmax=10
plt.figure(figsize=(15,5))
plt.subplot(121)
plt.plot(up2[:,0],up2[:,2],',',color='C1')
plt.plot(up1[:,0],up1[:,2],',',color='C0')
plt.fill([xmin,xmax,xmax,xmin],[-1,-1,1,1],color="k",alpha=0.4)
plt.ylim(-1,1)

ene1=1.0/np.sqrt(1.0-up1[:,2]**2-up1[:,3]**2-up1[:,4]**2)-1
ene2=1.0/np.sqrt(1.0-up2[:,2]**2-up2[:,3]**2-up2[:,4]**2)-1

plt.subplot(122)
plt.hist(ene1[np.where((up1[:,0]>=xmin)&(up1[:,0]<xmax))],bins=np.logspace(np.log10(enmin),np.log10(enmax),bins),log=True,histtype='step',color='C0')
#plt.hist(ene1[:],bins=np.logspace(np.log10(1e-6),np.log10(1.0),200),log=True,histtype='step',color='b')
plt.hist(ene2[np.where((up2[:,0]>=xmin)&(up2[:,0]<xmax))],bins=np.logspace(np.log10(enmin),np.log10(enmax),bins),log=True,histtype='step',color='C1')
#plt.hist(ene2[:],bins=np.logspace(np.log10(1e-6),np.log10(1.0),200),log=True,histtype='step',color='orange')
plt.yscale('log')
plt.xscale('log')
plt.xlim(enmin,enmax)
plt.show()

Animation

In [ ]:
%%time
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np

fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(15,5))
plt.subplots_adjust(hspace=0.4)

def update_anim(it):
    for ax in (ax1, ax2, ax3, ax4):
        ax.clear()
    iit=intvl*(it+1)
    nmv=['ex','ey','ez','bx','by','bz']
    EB=np.zeros((6,(ny+2)*nrank,nx+2))
    for iit in range(35,45):
        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+2)
                EB[i,(ny+2)*irank:(ny+2)*(irank+1),:]=d
            
    im1=ax1.imshow(EB[0,:,:],aspect='auto',origin='lower',cmap='terrain_r');ax1.set_title(r'$B_x$')
    im2=ax2.imshow(EB[1,:,:],aspect='auto',origin='lower',cmap='terrain_r');ax2.set_title(r'$B_y$')
    im3=ax3.imshow(EB[2,:,:],aspect='auto',origin='lower',cmap='terrain_r');ax3.set_title(r'$B_z$')
    im4=ax4.imshow(EB[3,:,:],aspect='auto',origin='lower',cmap='terrain_r');ax4.set_title(r'$E_x$')
    im5=ax5.imshow(EB[4,:,:],aspect='auto',origin='lower',cmap='terrain_r');ax5.set_title(r'$E_y$')
    im6=ax6.imshow(EB[5,:,:],aspect='auto',origin='lower',cmap='terrain_r');ax6.set_title(r'$E_z$')
plt.rcParams['animation.html'] = 'html5'
anim = FuncAnimation(fig, update_anim, blit=False,frames=41)
plt.close()
In [ ]:
%%time
anim

Animation: another method (output png files)

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+2)*nrank,nx+2))
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+2)
            EB[i,(ny+2)*irank:(ny+2)*(irank+1),:]=d

    fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(16,8))
    plt.subplots_adjust(hspace=0.5,wspace=0.3)
    fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt), fontsize=16)

    im1=ax1.imshow(EB[3,:,:],aspect='auto',origin='lower',cmap='terrain_r',clim=(-5,5));plt.colorbar(im1,ax=ax1);ax1.set_title(r'$B_x$')
    im2=ax2.imshow(EB[4,:,:],aspect='auto',origin='lower',cmap='terrain_r',clim=(0,10));plt.colorbar(im2,ax=ax2);ax2.set_title(r'$B_y$')
    im3=ax3.imshow(EB[5,:,:],aspect='auto',origin='lower',cmap='terrain_r',clim=(-5,5));plt.colorbar(im3,ax=ax3);ax3.set_title(r'$B_z$')
    im4=ax4.imshow(EB[0,:,:],aspect='auto',origin='lower',cmap='terrain_r',clim=(-1,1));plt.colorbar(im4,ax=ax4);ax4.set_title(r'$E_x$')
    im5=ax5.imshow(EB[1,:,:],aspect='auto',origin='lower',cmap='terrain_r',clim=(-1,1));plt.colorbar(im5,ax=ax5);ax5.set_title(r'$E_y$')
    im6=ax6.imshow(EB[2,:,:],aspect='auto',origin='lower',cmap='terrain_r',clim=(-.5,.5));plt.colorbar(im6,ax=ax6);ax6.set_title(r'$E_z$')
    
    #ax1.plot(np.mean(EB[3,:,:],axis=0)*5+160,color='k')
    ax2.plot(np.mean(EB[4,:,:],axis=0)*10,color='w')
    #ax3.plot(np.mean(EB[5,:,:],axis=0)*5+160,color='k')
    #ax4.plot(np.mean(EB[0,:,:],axis=0)*5+160,color='k')
    #ax5.plot(np.mean(EB[1,:,:],axis=0)*5+160,color='k')
    #ax6.plot(np.mean(EB[2,:,:],axis=0)*5+160,color='k')
    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$')

    fig.savefig('figures/EB%03d.png' %(iit), bbox_inches='tight')
    #print(iit)
    plt.close()
In [ ]:
!convert -delay 10 -loop 0 figures/EB*.png movie.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