Electrom Beam Instability

We show an example of the electron beam instability. The main calculation has been done by 1D full-PIC simulation written in FORTRAN.

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

module input
!### definition by a user ###!
integer, parameter :: ppc=200
integer, parameter :: nx=512
integer, parameter :: nsp=3
integer, parameter :: nt=4096
integer, parameter :: isavw=2, isavp=128
real(8), parameter :: theta=0.505d0
real(8), parameter :: c =1.d0
real(8), parameter :: dx=0.05d0
real(8), parameter :: dt=0.05d0
real(8), parameter :: b0=0.1d0
real(8), parameter :: th=0
real(8), parameter :: q(nsp)    =(/ -1.d0             , 1.d0             , -1.d0              /)
real(8), parameter :: rm(nsp)   =(/  1.d0             , 100.d0            , 1.d0             /)
real(8), parameter :: qomdt(nsp)=(/  0.5*dt*q(1)/rm(1), 0.5*dt*q(2)/rm(2), 0.5*dt*q(3)/rm(3) /)
real(8), parameter :: dn(nsp)   =(/  0.9d0/ppc         , 1.d0/ppc         , 0.1d0/ppc          /)
real(8), parameter :: vt(nsp)   =(/  0.05d0           , 0.05d0/sqrt(rm(2))     , 0.01d0/rm(3)        /)
real(8), parameter :: vx0(nsp)  =(/  -0.02d0          , 0.00d0           , 0.18d0              /)

!### don't change if you have no reason! ###!
integer, parameter :: npmax=ppc*nx
integer, dimension(nsp) :: np
integer, dimension(nx) :: iip,iim
real(8), parameter :: pi=3.141592653589793d0
real(8), parameter :: lx=dx*nx



end module input
In [2]:
%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 [3]:
nx=512
nd=2048   #the number of output files for EM fields
isav=2
dx=0.05
dt=0.05*isav
lx=nx*dx
nb=0.1

Space-time evolution of electromagnetic fields

In [4]:
f1 =open('data/exhst.d'  ,'rb')
f2 =open('data/eyhst.d'  ,'rb')
f3 =open('data/ezhst.d'  ,'rb')
f4 =open('data/bxhst.d'  ,'rb')
f5 =open('data/byhst.d'  ,'rb')
f6 =open('data/bzhst.d'  ,'rb')
f7 =open('data/rhohst.d' ,'rb')
f8 =open('data/jxhst.d'  ,'rb')
f9 =open('data/jyhst.d'  ,'rb')
f10=open('data/jzhst.d'  ,'rb')

exhst =np.fromfile(f1, dtype='float64').reshape(-1,nx)
eyhst =np.fromfile(f2, dtype='float64').reshape(-1,nx)
ezhst =np.fromfile(f3, dtype='float64').reshape(-1,nx)
bxhst =np.fromfile(f4, dtype='float64').reshape(-1,nx)
byhst =np.fromfile(f5, dtype='float64').reshape(-1,nx)
bzhst =np.fromfile(f6, dtype='float64').reshape(-1,nx)
rhohst=np.fromfile(f7, dtype='float64').reshape(-1,nx)
jxhst =np.fromfile(f8, dtype='float64').reshape(-1,nx)
jyhst =np.fromfile(f9, dtype='float64').reshape(-1,nx)
jzhst =np.fromfile(f10,dtype='float64').reshape(-1,nx)
In [5]:
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(235, sharex=ax1, sharey=ax1)
ax5 = fig.add_subplot(236, sharex=ax1, sharey=ax1)
im1=ax1.imshow(exhst,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,lx,0,tmax])
im2=ax2.imshow(eyhst,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,tmax],clim=(-0.001,0.001))
im3=ax3.imshow(ezhst,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,tmax],clim=(-0.001,0.001))
im4=ax4.imshow(byhst,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,tmax],clim=(-0.001,0.001))
im5=ax5.imshow(bzhst,aspect='auto',origin='lower',cmap='bwr'      ,extent=[0,lx,0,tmax],clim=(-0.001,0.001))
ax1.set_xlabel(r'$x/(c/\omega_{e})$')
ax4.set_xlabel(r'$x/(c/\omega_{e})$')
ax5.set_xlabel(r'$x/(c/\omega_{e})$')
ax1.set_ylabel(r'$t\omega_e$')
ax4.set_ylabel(r'$t\omega_e$')
ax1.set_xlim(0,lx);
ax1.set_title(r'$E_x$')
ax2.set_title(r'$E_y$')
ax3.set_title(r'$E_z$')
ax4.set_title(r'$B_y$')
ax5.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.show()

Time evolution of electromagnetic energy

In [6]:
tt=dt*np.arange(nd)
enex=np.zeros((nd))
eney=np.zeros((nd))
enez=np.zeros((nd))
enby=np.zeros((nd))
enbz=np.zeros((nd))
for it in range(nd):
    enex[it]=np.average(exhst[it,:]**2)
    eney[it]=np.average(eyhst[it,:]**2)
    enez[it]=np.average(ezhst[it,:]**2)
    enby[it]=np.average(byhst[it,:]**2)
    enbz[it]=np.average(bzhst[it,:]**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()
plt.xlabel(r'$t\omega_e$')
plt.yscale('log')
plt.ylim(1e-10,1e-2)
plt.show()

Dispersion Relation

Note: FFT2 of Scipy is definec by
$$y(j) = \sum\limits_{k = 0 \ldots n - 1} {x[k]*\exp \left( { - \sqrt { - 1} *j*k*2*\pi /n} \right)} ,\,\,j = 0 \ldots n - 1$$ for the both axes.

Therefore, if you would like to calculate waves of $\exp \left[ {\sqrt { - 1} \left( {kx - \omega t} \right)} \right]$, you have to swritch the column or raw of a matrix.

In [7]:
kmin=2*mt.pi/(dx*nx)*(-nx/2); kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nd)*(-nd/2); wmax=2*mt.pi/(dt*nd)*(nd/2)
kaxis=np.linspace(kmin,kmax,nx); waxis=np.linspace(wmin,wmax,nd)

fftexsave=fftshift(fft2(exhst[:,::-1]))  #switch the raw
plt.imshow(np.abs(fftexsave),interpolation='none',origin='lower',aspect='auto',cmap='terrain',extent=([kmin,kmax,wmin,wmax]))    
plt.axis([kmin/2,kmax/2,0,wmax/2])
plt.colorbar()
plt.clim(0,150)
plt.xlabel(r'$k(c/\omega_{pe})$')
plt.ylabel(r'$\omega/\omega_{pe}$')
plt.show()

Animation

In [9]:
bs=256

plt.rcParams['animation.html'] = 'html5'
#fig,(ax1,ax2)=plt.subplots(1,3,figsize=(10,5))
fig=plt.figure(figsize=(12,4))
ax1=plt.subplot2grid((1,3),(0,0),colspan=2)
ax2=plt.subplot2grid((1,3),(0,2))

def update_anim(i):
    it=i+1
    f1=open('data/up1%05d.d' %(it),'rb') 
    f2=open('data/up2%05d.d' %(it),'rb')
    f3=open('data/up3%05d.d' %(it),'rb')
    xp1=np.fromfile(f1,dtype='float64').reshape(-1,5)
    xp2=np.fromfile(f2,dtype='float64').reshape(-1,5)
    xp3=np.fromfile(f3,dtype='float64').reshape(-1,5)
    
    for ax in (ax1,ax2): 
        ax.clear()

    ax1.plot(xp1[:,0],xp1[:,1],'C0,')
    ax1.plot(xp2[:,0],xp2[:,1],'C1,')
    ax1.plot(xp3[:,0],xp3[:,1],'C2,')
    ax1.set_xlabel(r'$x/(c/\omega_{pe})$')
    ax1.set_ylabel(r'$v_x/c$')
    ax1.set_ylim(-0.5,0.5)

    ax2.hist(xp1[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='b',weights=(1-nb)*np.ones(len(xp1)))
    ax2.hist(xp2[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='orange')
    ax2.hist(xp3[:,1],bins=bs,range=(-0.25,0.25),histtype='step',color='g',weights=nb*np.ones(len(xp1)))
    ax2.set_xlabel(r'$v_x/c$')
    ax2.set_ylabel(r'$f(v_x)$')
    #ax2.set_ylim(0,10000)
    
anim=animation.FuncAnimation(fig,update_anim,frames=32)    
plt.close()
anim
Out[9]: