Weibel Instability

We show an example of the electron beam instability. The main calculation has been done by 2D electron-PIC simulation written in FORTRAN. The electron-PIC simulation means that only electrons are solved and ions are fixed in a simulation box.

In [15]:
%%bash
cat input.f90

module input
!### definition by a user ###!
integer, parameter :: ppc=20
integer, parameter :: nx=256
integer, parameter :: ny=256
integer, parameter :: nsp=2
integer, parameter :: nt=4096
integer, parameter :: isavw=16, isavp=256
real(8), parameter :: theta=0.505d0
real(8), parameter :: c =1.d0
real(8), parameter :: dx=0.1d0
real(8), parameter :: dy=0.1d0
real(8), parameter :: dt=0.1d0
real(8), parameter :: b0=0.1d0
real(8), parameter :: th=0, phi=90
real(8), parameter :: q(nsp)    =(/ -1.d0             ,-1.d0  /)
real(8), parameter :: rm(nsp)   =(/  1.d0             , 1.d0  /)
real(8), parameter :: qomdt(nsp)=(/  0.5*dt*q(1)/rm(1), 0.5*dt*q(2)/rm(2)  /)
real(8), parameter :: dn(nsp)   =(/  0.5d0/ppc        , 0.5d0/ppc          /)
real(8), parameter :: vt(nsp)  =(/    0.1d0          ,  0.1d0/)
real(8), parameter :: vz0(nsp)  =(/   -0.3d0          ,  0.3d0/)

!### don't change if you have no reason! ###!
integer, parameter :: npmax=ppc*nx*ny
integer, dimension(nsp) :: np
integer, dimension(nx) :: iipx,iimx
integer, dimension(ny) :: iipy,iimy
integer, dimension(nx*ny) :: ijmx, ijpx, ijmy, ijpy
real(8), parameter :: pi=3.141592653589793d0
real(8), parameter :: lx=dx*nx, ly=dy*ny

end module input
In [16]:
%config InlineBackend.figure_format = 'retina'
from numba import jit,f8,u8
import numpy as np
import math as mt
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,fftshift
import time
from IPython import display
import matplotlib.animation as animation
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [19]:
nx=256
ny=256
nd=256
isav=16
dx=0.1; dy=0.1
dt=0.1*isav
lx=nx*dx; ly=ny*dy

Snapshot of electromagnetic fields

In [24]:
it=100
f1 =open('data/ex%05d.d' %(it) ,'rb')
f2 =open('data/ey%05d.d' %(it) ,'rb')
f3 =open('data/ez%05d.d' %(it) ,'rb')
f4 =open('data/bx%05d.d' %(it) ,'rb')
f5 =open('data/by%05d.d' %(it) ,'rb')
f6 =open('data/bz%05d.d' %(it) ,'rb')
f7 =open('data/rho%05d.d'%(it) ,'rb')
f8 =open('data/jx%05d.d' %(it) ,'rb')
f9 =open('data/jy%05d.d' %(it) ,'rb')
f10=open('data/jz%05d.d' %(it) ,'rb')

ex =np.fromfile(f1, dtype='float64' ).reshape(-1,nx)
ey =np.fromfile(f2, dtype='float64' ).reshape(-1,nx)
ez =np.fromfile(f3, dtype='float64' ).reshape(-1,nx)
bx =np.fromfile(f4, dtype='float64' ).reshape(-1,nx)
by =np.fromfile(f5, dtype='float64' ).reshape(-1,nx)
bz =np.fromfile(f6, dtype='float64' ).reshape(-1,nx)
rho=np.fromfile(f7, dtype='float64' ).reshape(-1,nx)
jx =np.fromfile(f8, dtype='float64' ).reshape(-1,nx)
jy =np.fromfile(f9, dtype='float64' ).reshape(-1,nx)
jz =np.fromfile(f10,dtype='float64' ).reshape(-1,nx)
tmax=nd*dt
fig=plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(231)
ax2 = fig.add_subplot(232, sharex=ax1, sharey=ax1)
ax3 = fig.add_subplot(233, sharex=ax1, sharey=ax1)
ax4 = fig.add_subplot(234, sharex=ax1, sharey=ax1)
ax5 = fig.add_subplot(235, sharex=ax1, sharey=ax1)
ax6 = fig.add_subplot(236, sharex=ax1, sharey=ax1)
im1=ax1.imshow(ex,aspect='auto',origin='lower',cmap='bwr',extent=[0,lx,0,ly])
im2=ax2.imshow(ey,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,ly])
im3=ax3.imshow(ez,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,ly])
im4=ax4.imshow(bx,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,ly])
im5=ax5.imshow(by,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,ly])
im6=ax6.imshow(bz,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,ly])
ax4.set_xlabel(r'$x/(c/\omega_{pe})$')
ax5.set_xlabel(r'$x/(c/\omega_{pe})$')
ax6.set_xlabel(r'$x/(c/\omega_{pe})$')
ax1.set_ylabel(r'$y/(c/\omega_{pe})$')
ax4.set_ylabel(r'$y/(c/\omega_{pe})$')
ax1.set_title(r'$E_x$')
ax2.set_title(r'$E_y$')
ax3.set_title(r'$E_z$')
ax4.set_title(r'$B_x$')
ax5.set_title(r'$B_y$')
ax6.set_title(r'$B_z$')
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)
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt))
plt.show()

Time evolution of electromagnetic fields

In [39]:
plt.figure(figsize=(8,4))
tt=dt*np.arange(nd)
enex=np.zeros((nd))
eney=np.zeros((nd))
enez=np.zeros((nd))
enbx=np.zeros((nd))
enby=np.zeros((nd))
enbz=np.zeros((nd))
for it in range(nd):
    f1 =open('data/ex%05d.d' %(it) ,'rb')
    f2 =open('data/ey%05d.d' %(it) ,'rb')
    f3 =open('data/ez%05d.d' %(it) ,'rb')
    f4 =open('data/bx%05d.d' %(it) ,'rb')
    f5 =open('data/by%05d.d' %(it) ,'rb')
    f6 =open('data/bz%05d.d' %(it) ,'rb')
    
    ex =np.fromfile(f1, dtype='float64' ).reshape(-1,nx)
    ey =np.fromfile(f2, dtype='float64' ).reshape(-1,nx)
    ez =np.fromfile(f3, dtype='float64' ).reshape(-1,nx)
    bx =np.fromfile(f4, dtype='float64' ).reshape(-1,nx)
    by =np.fromfile(f5, dtype='float64' ).reshape(-1,nx)
    bz =np.fromfile(f6, dtype='float64' ).reshape(-1,nx)
    enex[it]=np.average(ex**2)
    eney[it]=np.average(ey**2)
    enez[it]=np.average(ez**2)
    enbx[it]=np.average(bx**2)
    enby[it]=np.average(by**2)
    enbz[it]=np.average(bz**2)

plt.plot(tt,enex,label=r'$E_x^2$')
plt.plot(tt,eney,label=r'$E_y^2$')
plt.plot(tt,enez,label=r'$E_z^2$')
plt.plot(tt,enby,label=r'$B_y^2$')
plt.plot(tt,enbz,label=r'$B_z^2$')
plt.legend(loc='lower right')
plt.xlabel(r'$t\omega_e$')
plt.yscale('log')
plt.ylim(1e-6,1e-1)
plt.show()

Animation

In [41]:
plt.rcParams['animation.html'] = 'html5'
#fig,(ax1,ax2)=plt.subplots(1,3,figsize=(10,5))
fig=plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)

def update_anim(i):
    it=i*2
    #f1=open('data/rho%05d.d' %(it),'rb') 
    f2=open('data/bx%05d.d' %(it),'rb') 
    f3=open('data/by%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)
    
    #for ax in (ax1,ax2): 
    ax.clear()

    ax.imshow(np.sqrt(bx**2+by**2),origin='lower',aspect='auto',cmap='gist_heat')
    #ax.imshow(rho,origin='lower',aspect='auto',cmap='gist_heat')
    ax.set_title(r"$t\omega_{pe}=%.2f$" %(it*dt))
    ax.set_xlabel(r'$x/(c/\omega_{pe})$')
    ax.set_ylabel(r'$y/(c/\omega_{pe})$')
    
anim=animation.FuncAnimation(fig,update_anim,frames=128)    
plt.close()
anim
Out[41]: