Electron Beam 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 [1]:
%%bash
cat input.f90
In [5]:
%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 [10]:
nx=256
ny=128
nd=512
isav=8
dx=0.05; dy=0.05
dt=0.05*isav
lx=nx*dx; ly=ny*dy
nb=0.1
Space-time evolution of electromagnetic fields¶
Note: Integration for the $y$-axis.
In [11]:
exhst=np.zeros((nd,nx));eyhst=np.zeros((nd,nx));ezhst=np.zeros((nd,nx))
bxhst=np.zeros((nd,nx));byhst=np.zeros((nd,nx));bzhst=np.zeros((nd,nx))
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')
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)
exhst[it,:]=np.mean(ex,axis=0);eyhst[it,:]=np.mean(ey,axis=0);ezhst[it,:]=np.mean(ez,axis=0)
bxhst[it,:]=np.mean(bx,axis=0);byhst[it,:]=np.mean(by,axis=0);bzhst[it,:]=np.mean(bz,axis=0)
In [12]:
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(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])
im3=ax3.imshow(ezhst,aspect='auto',origin='lower',cmap='bwr' ,extent=[0,lx,0,tmax])
im4=ax4.imshow(bxhst,aspect='auto',origin='lower',cmap='bwr' ,extent=[0,lx,0,tmax])
im5=ax5.imshow(byhst,aspect='auto',origin='lower',cmap='bwr' ,extent=[0,lx,0,tmax])
im6=ax6.imshow(bzhst,aspect='auto',origin='lower',cmap='bwr' ,extent=[0,lx,0,tmax])
ax4.set_xlabel(r'$x/(c/\omega_{e})$')
ax5.set_xlabel(r'$x/(c/\omega_{e})$')
ax6.set_xlabel(r'$x/(c/\omega_{e})$')
ax1.set_ylabel(r'$t\omega_e$')
ax4.set_ylabel(r'$t\omega_e$')
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(im5,ax=ax6)
plt.show()
Time evolution of electromagnetic fields¶
In [13]:
plt.figure(figsize=(8,4))
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(loc='lower right')
plt.xlabel(r'$t\omega_e$')
plt.yscale('log')
plt.ylim(1e-10,1e-2)
plt.show()
Dispersion relation¶
In [14]:
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)
fftex=fftshift(fft2(exhst[:,::-1]))
fftby=fftshift(fft2(byhst[:,::-1]))
plt.figure(figsize=(12,12))
plt.subplot(221);plt.imshow(exhst,origin='lower',cmap='terrain_r',aspect='auto',extent=(0,nx*dx,0,nd*dt))
plt.title(r'$E_x$');plt.xlabel(r'$x/(c/\omega_{pe})$');plt.ylabel(r'$t\omega_{pe}$')
plt.subplot(222);plt.imshow(byhst,origin='lower',cmap='bwr',aspect='auto',extent=(0,nx*dx,0,nd*dt))
plt.title(r'$B_z$');plt.xlabel(r'$x/(c/\omega_{pe})$');plt.ylabel(r'$t\omega_{pe}$')
plt.subplot(223);plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.axis([kmin,kmax,0,wmax/2])
plt.xlabel(r'$k(c/\omega_{pe})$');plt.ylabel(r'$\omega/\omega_{pe}$')
plt.clim(0,1)
plt.subplot(224);plt.imshow(abs(fftby),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.xlabel(r'$k(c/\omega_{pe})$');plt.ylabel(r'$\omega/\omega_{pe}$')
plt.axis([kmin/8,kmax/8,0,wmax/2])
plt.clim(0,1)
plt.show()
Snapshot of electromagnetic fields¶
In [15]:
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='terrain_r',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(im5,ax=ax6)
fig.suptitle(r"$t\omega_{pe}=%.2f$" %(it*dt))
plt.show()
Animation¶
In [16]:
plt.rcParams['animation.html'] = 'html5'
#fig,(ax1,ax2)=plt.subplots(1,3,figsize=(10,5))
fig=plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
def update_anim(i):
it=i*4
f1=open('data/rho%05d.d' %(it),'rb')
rho=np.fromfile(f1,dtype='float64').reshape(-1,nx)
#for ax in (ax1,ax2):
ax.clear()
ax.imshow(rho,origin='lower',aspect='auto')
ax.set_title(r"$t\omega_{pe}=%.2f$" %(it*dt))
ax.set_xlabel(r'$x/(c/\omega_{pe})$')
ax.set_ylabel(r'$v_x/c$')
anim=animation.FuncAnimation(fig,update_anim,frames=128)
plt.close()
anim
Out[16]: