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 the simulation, ions are fixed in the simulation box.
In [1]:
%%bash
cat input.f90
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 [13]:
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')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,5)
xp2=np.fromfile(f2,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.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',weights=nb*np.ones(len(xp2)))
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[13]: