Linear Waves

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

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

module input
!### definition by a user ###!
integer, parameter :: ppc=20
integer, parameter :: nx=512
integer, parameter :: ny=512
integer, parameter :: nsp=2
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 :: dy=0.05d0
real(8), parameter :: dt=0.05d0
real(8), parameter :: b0=2.d0
real(8), parameter :: th=0, phi=0
real(8), parameter :: q(nsp)    =(/ -1.d0             , 1.d0              /)
real(8), parameter :: rm(nsp)   =(/  1.d0             , 4.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)   =(/  1.d0/ppc         , 1.d0/ppc          /)
real(8), parameter :: vt(nsp)  =(/  0.05d0             , 0.05d0/sqrt(rm(2)) /)
real(8), parameter :: vx0(nsp)  =(/  0.d0             , 0.d0              /)

!### 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 [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
ny=512
nd=2048
isav=2
dx=0.05; dy=0.05
dt=0.05*isav
lx=nx*dx; ly=ny*dy
lt=dt*nd
In [4]:
exhstx=np.zeros((nd,nx));eyhstx=np.zeros((nd,nx));ezhstx=np.zeros((nd,nx))
bxhstx=np.zeros((nd,nx));byhstx=np.zeros((nd,nx));bzhstx=np.zeros((nd,nx))
exhsty=np.zeros((nd,ny));eyhsty=np.zeros((nd,ny));ezhsty=np.zeros((nd,ny))
bxhsty=np.zeros((nd,ny));byhsty=np.zeros((nd,ny));bzhsty=np.zeros((nd,ny))
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)
    
    exhstx[it,:]=np.mean(ex,axis=0);eyhstx[it,:]=np.mean(ey,axis=0);ezhstx[it,:]=np.mean(ez,axis=0)
    bxhstx[it,:]=np.mean(bx,axis=0);byhstx[it,:]=np.mean(by,axis=0);bzhstx[it,:]=np.mean(bz,axis=0)
    exhsty[it,:]=np.mean(ex,axis=1);eyhsty[it,:]=np.mean(ey,axis=1);ezhsty[it,:]=np.mean(ez,axis=1)
    bxhsty[it,:]=np.mean(bx,axis=1);byhsty[it,:]=np.mean(by,axis=1);bzhsty[it,:]=np.mean(bz,axis=1)

Dispersion relation along $x-$axis ($\theta=0^\circ$)

In [5]:
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(exhstx[:,::-1]))
fftey=fftshift(fft2(eyhstx[:,::-1]))
fftez=fftshift(fft2(ezhstx[:,::-1]))
fftbx=fftshift(fft2(bxhstx[:,::-1]))
fftby=fftshift(fft2(byhstx[:,::-1]))
fftbz=fftshift(fft2(bzhstx[:,::-1]))

fig, ((ax1,ax2,ax3,ax4),(ax5,ax6,ax7,ax8),(ax9,ax10,ax11,ax12))=plt.subplots(3,4,figsize=(25,15))
ax1. imshow(exhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax5. imshow(eyhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax9. imshow(ezhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax3. imshow(bxhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax7. imshow(byhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax11.imshow(bzhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax1 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax5 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax9 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax3 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax7 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax11.set_xlabel(r'$x/(c/\omega_{pe})$');
ax1 .set_ylabel(r'$t\omega_{pe}$')
ax5 .set_ylabel(r'$t\omega_{pe}$')
ax9 .set_ylabel(r'$t\omega_{pe}$')
ax3 .set_ylabel(r'$t\omega_{pe}$')
ax7 .set_ylabel(r'$t\omega_{pe}$')
ax11.set_ylabel(r'$t\omega_{pe}$')

ax2. imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax6. imshow(abs(fftey),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax10.imshow(abs(fftez),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax4. imshow(abs(fftbx),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax8. imshow(abs(fftby),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax12.imshow(abs(fftbz),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax2 .set_xlabel(r'$kc/\omega_{pe}$');
ax6 .set_xlabel(r'$kc/\omega_{pe}$');
ax10.set_xlabel(r'$kc/\omega_{pe}$');
ax4 .set_xlabel(r'$kc/\omega_{pe}$');
ax8 .set_xlabel(r'$kc/\omega_{pe}$');
ax12.set_xlabel(r'$kc/\omega_{pe}$');
ax2 .set_ylabel(r'$\omega/\omega_{pe}$')
ax6 .set_ylabel(r'$\omega/\omega_{pe}$')
ax10.set_ylabel(r'$\omega/\omega_{pe}$')
ax4 .set_ylabel(r'$\omega/\omega_{pe}$')
ax8 .set_ylabel(r'$\omega/\omega_{pe}$')
ax12.set_ylabel(r'$\omega/\omega_{pe}$')

ax2 .axis([0,kmax/4,0,wmax/4])
ax6 .axis([0,kmax/4,0,wmax/4])
ax10.axis([0,kmax/4,0,wmax/4])
ax4 .axis([0,kmax/4,0,wmax/4])
ax8 .axis([0,kmax/4,0,wmax/4])
ax12.axis([0,kmax/4,0,wmax/4])

plt.show()

Dispersion relation along $y-$axis ($\theta=90^\circ$)

In [7]:
kmin=2*mt.pi/(dy*ny)*(-ny/2); kmax=2*mt.pi/(dy*ny)*(ny/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(exhsty[:,::-1]))
fftey=fftshift(fft2(eyhsty[:,::-1]))
fftez=fftshift(fft2(ezhsty[:,::-1]))
fftbx=fftshift(fft2(bxhsty[:,::-1]))
fftby=fftshift(fft2(byhsty[:,::-1]))
fftbz=fftshift(fft2(bzhsty[:,::-1]))

fig, ((ax1,ax2,ax3,ax4),(ax5,ax6,ax7,ax8),(ax9,ax10,ax11,ax12))=plt.subplots(3,4,figsize=(25,15))
ax1. imshow(exhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax5. imshow(eyhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax9. imshow(ezhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax3. imshow(bxhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax7. imshow(byhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax11.imshow(bzhstx,origin='lower',cmap='bwr',aspect='auto',extent=[0,lx,0,lt]);
ax1 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax5 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax9 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax3 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax7 .set_xlabel(r'$x/(c/\omega_{pe})$');
ax11.set_xlabel(r'$x/(c/\omega_{pe})$');
ax1 .set_ylabel(r'$t\omega_{pe}$')
ax5 .set_ylabel(r'$t\omega_{pe}$')
ax9 .set_ylabel(r'$t\omega_{pe}$')
ax3 .set_ylabel(r'$t\omega_{pe}$')
ax7 .set_ylabel(r'$t\omega_{pe}$')
ax11.set_ylabel(r'$t\omega_{pe}$')

ax2. imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax6. imshow(abs(fftey),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax10.imshow(abs(fftez),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax4. imshow(abs(fftbx),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax8. imshow(abs(fftby),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax12.imshow(abs(fftbz),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]),clim=(0,1));
ax2 .set_xlabel(r'$kc/\omega_{pe}$');
ax6 .set_xlabel(r'$kc/\omega_{pe}$');
ax10.set_xlabel(r'$kc/\omega_{pe}$');
ax4 .set_xlabel(r'$kc/\omega_{pe}$');
ax8 .set_xlabel(r'$kc/\omega_{pe}$');
ax12.set_xlabel(r'$kc/\omega_{pe}$');
ax2 .set_ylabel(r'$\omega/\omega_{pe}$')
ax6 .set_ylabel(r'$\omega/\omega_{pe}$')
ax10.set_ylabel(r'$\omega/\omega_{pe}$')
ax4 .set_ylabel(r'$\omega/\omega_{pe}$')
ax8 .set_ylabel(r'$\omega/\omega_{pe}$')
ax12.set_ylabel(r'$\omega/\omega_{pe}$')

ax2 .axis([0,kmax/4,0,wmax/4])
ax6 .axis([0,kmax/4,0,wmax/4])
ax10.axis([0,kmax/4,0,wmax/4])
ax4 .axis([0,kmax/4,0,wmax/4])
ax8 .axis([0,kmax/4,0,wmax/4])
ax12.axis([0,kmax/4,0,wmax/4])

plt.show()