In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from numba.decorators import jit
from numba import f8,u8
import scipy.fftpack as sf
from IPython.display import display, clear_output

Electrostatic 1-D Particle in Cell simulation

In [2]:
def es1(vtc,vtb,ppc,dx,nx,dt,nt,nb,vb):
    np.seterr(divide='ignore', invalid='ignore')
    nop=ppc*nx
    lx=dx*nx
    xx=dx*np.arange(nx)
    tt=dt*np.arange(nt)
    xp=np.zeros((2,nop))
    vp=np.zeros((2,nop))
    xp[0,:]=np.linspace(0,lx-lx/nop,nop)
    xp[1,:]=np.linspace(0,lx-lx/nop,nop)
    vp[0,:]=np.random.normal(0,vtc,nop)-vb*nb/(1-nb)
    vp[1,:]=np.random.normal(0,vtb,nop)+vb

    ds=np.zeros((2,nx))
    kk=2*mt.pi/lx*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    Kl=(np.sin(0.5*kk*dx)/(0.5*dx))**2
    kl=np.sin(kk*dx)/dx
    
    esave=np.zeros((nt,nx))
    
    print('initial condition')
    plt.subplot(121);plt.plot(xp[0,:],vp[0,:],',');plt.plot(xp[1,:],vp[1,:],',')
    plt.subplot(122);plt.hist(vp[0,:],bins=500,histtype='step',weights=(1-nb)*np.ones(nop)); plt.hist(vp[1,:], bins=500, histtype='step', weights=nb*np.ones(nop))
    plt.show()
    
    q   =np.zeros(2) 
    q[0]=nx/nop*(1-nb)
    q[1]=nx/nop*nb    
    
    for it in range(nt):
        xp=xp+dt*vp
        xp[xp>lx]=xp[xp>lx]-lx
        xp[xp<0 ]=xp[xp<0 ]+lx
        
        ds=np.zeros((2,nx))
        ds=dens(dx,nx,nop,q,xp,ds)
        #exfft=1j/kk*sf.fft(np.sum(ds,0))
        #exfft[0]=0
        exfft=1j*kl*sf.fft(np.sum(ds,0))/Kl
        exfft[0]=0
        ex=np.real(sf.ifft(exfft))        
        
        vp=acc(dx,nx,nop,xp,vp,ex,dt)
        esave[it,:]=ex[:]
        
    kmin=2*mt.pi/(dx*nx)*(-nx/2); kmax=2*mt.pi/(dx*nx)*(nx/2)
    wmin=2*mt.pi/(dt*nt)*(-nt/2); wmax=2*mt.pi/(dt*nt)*(nt/2)
    kaxis=np.linspace(kmin,kmax,nx); waxis=np.linspace(wmin,wmax,nt)
    
    return locals()

@jit('f8[:,:](f8,u8,u8,f8[:,:],f8[:,:],f8[:],f8)')
def acc(dx,nx,nop,xp,vp,ex,dt):
    for isp in range(2):
        for ip in range(nop):
            ixm=mt.floor(xp[isp,ip]/dx); ixp=ixm+1
            wxp=xp[isp,ip]/dx-ixm; wxm=1-wxp

            if ixp>nx-1: ixp=ixp-nx
            vp[isp,ip]=vp[isp,ip]-dt*(wxm*ex[ixm]+wxp*ex[ixp])
        
    return vp

@jit('f8[:,:](f8,u8,u8,f8[:],f8[:,:],f8[:,:])')
def dens(dx,nx,nop,q,xp,ds):
    for isp in range(2):
        for ip in range(nop):
            ixm=mt.floor(xp[isp,ip]/dx); ixp=ixm+1
            wxp=xp[isp,ip]/dx-ixm; wxm=1-wxp

            if ixp>nx-1: ixp=ixp-nx

            ds[isp,ixm]=ds[isp,ixm]+wxm
            ds[isp,ixp]=ds[isp,ixp]+wxp 
        ds[isp,:]=ds[isp,:]*q[isp]
    return ds

Linear waves

note: If you do not consider a beam ($v_b=0$), you have to set $n_b$ finite.

In [3]:
%%time 
data=es1(0.05,0.05,400,0.1,256,0.1,512,0.5,0) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 5.16 s, sys: 203 ms, total: 5.36 s
Wall time: 7.12 s
Space time evolution of $E_x$ and the dispersion relation
In [4]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
plt.clim(0,10)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,1)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],wmin,wmax]))
plt.axis([0,tt[nt-1],wmin/4,wmax/4])
plt.clim(0,0.1)

plt.show()

Two fluid instability $n_b=0.01$, $v_b/v_{thc}=2$

In [5]:
%time data=es1(0.05,0.01,400,0.2,2**8,0.2,1024,0.01, 0.1) #(vtc,vtb,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.05 s, sys: 109 ms, total: 9.16 s
Wall time: 11 s
Space time evolution of $E_x$ and the dispersion relation
In [6]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.plot(kaxis,vb*kaxis,'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
plt.clim(0,10)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,1)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],wmin,wmax]))
plt.axis([0,tt[nt-1],wmin/4,wmax/4])
plt.clim(0,.1)
plt.show()

Two fluid instability $n_b=0.01$, $v_b/v_{thc}=5$

In [7]:
%time data=es1(0.05,0.01,400,0.2,2**8,0.2,1024,0.01, 0.25) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.53 s, sys: 62.5 ms, total: 9.59 s
Wall time: 13.6 s
Space time evolution of $E_x$ and the dispersion relation
In [8]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.plot(kaxis,vb*kaxis,'w--')
plt.axis([kmin/2,kmax/2,wmin/2,wmax/2])
plt.clim(0,50)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,1)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],wmin,wmax]))
plt.axis([0,tt[nt-1],wmin/2,wmax/2])

plt.show()

Two fluid instability $n_b=0.01$, $v_b/v_{thc}=10$

In [9]:
%time data=es1(0.05,0.01,400,0.2,2**8,0.2,1024,0.01, 0.5) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.61 s, sys: 219 ms, total: 9.83 s
Wall time: 14.9 s
Space time evolution of $E_x$ and the dispersion relation
In [10]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.plot(kaxis,vb*kaxis,'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
plt.clim(0,100)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,10)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],kmin,kmax]))
plt.axis([0,tt[nt-1],kmin/2,kmax/2])
plt.clim(0,1)
plt.show()

Two fluid instability $n_b=0.1$, $v_b/v_{thc}=2$

In [11]:
%time data=es1(0.05,0.01,400,0.2,2**8,0.2,1024,0.1, 0.1) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 8.58 s, sys: 234 ms, total: 8.81 s
Wall time: 12.1 s
Space time evolution of $E_x$ and the dispersion relation
In [12]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.plot(kaxis,vb*kaxis,'w--')
plt.axis([kmin,kmax,wmin/4,wmax/4])
plt.clim(0,10)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,1)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],kmin,kmax]))
plt.axis([0,tt[nt-1],kmin/2,kmax/2])
plt.clim(0,.1)
plt.show()

Two fluid instability $n_b=0.1$, $v_b/v_{thc}=5$

In [13]:
%time data=es1(0.05,0.01,400,0.2,2**8,0.2,1024,0.1, 0.25) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.08 s, sys: 188 ms, total: 9.27 s
Wall time: 11.4 s
Space time evolution of $E_x$ and the dispersion relation
In [14]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.plot(kaxis,vb*kaxis,'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
plt.clim(0,100)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,1)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],kmin,kmax]))
plt.axis([0,tt[nt-1],kmin,kmax])
plt.clim(0,0.5)
plt.show()

Two fluid instability $n_b=0.1$, $v_b/v_{thc}=10$

In [15]:
%time data=es1(0.05,0.01,400,0.2,2**8,0.2,1024,0.1, 0.5) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.75 s, sys: 156 ms, total: 9.91 s
Wall time: 13 s
Space time evolution of $E_x$ and the dispersion relation
In [16]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.plot(kaxis,vb*kaxis,'w--')
plt.axis([kmin/2,kmax/2,wmin/4,wmax/4])
plt.clim(0,50)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
plt.clim(0,1)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],kmin,kmax]))
plt.axis([0,tt[nt-1],kmin,kmax])
plt.clim(0,0.5)
plt.show()

Two fluid instability $n_b=0.5$, $v_b/v_{thc}=2$

In [17]:
%time data=es1(0.05,0.05,400,0.2,2**8,0.2,1024,0.5, 0.1) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.48 s, sys: 62.5 ms, total: 9.55 s
Wall time: 11.8 s
Space time evolution of $E_x$ and the dispersion relation
In [18]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
#plt.clim(0,50)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/10,wmax/10])
#plt.clim(0,10)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],wmin,wmax]))
plt.axis([0,tt[nt-1],wmin/2,wmax/2])

plt.show()

Two fluid instability $n_b=0.5$, $v_b/v_{thc}=5$

In [19]:
%time data=es1(0.05,0.05,400,0.2,2**8,0.2,1024,0.5, 0.25) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 9.72 s, sys: 109 ms, total: 9.83 s
Wall time: 13.3 s
Space time evolution of $E_x$ and the dispersion relation
In [20]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
#plt.clim(0,500)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
#plt.clim(0,10)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],wmin,wmax]))
plt.axis([0,tt[nt-1],wmin/4,wmax/4])

plt.show()

Two fluid instability $n_b=0.5$, $v_b/v_{thc}=10$

In [21]:
%%time
data=es1(0.05,0.05,400,0.2,2**8,0.2,1024,0.5, 0.5) #(vt,nop,dx,nx,dt,nt,nb,vb)
locals().update(data)

initial condition

CPU times: user 10.3 s, sys: 203 ms, total: 10.5 s
Wall time: 17.2 s
Space time evolution of $E_x$ and the dispersion relation
In [22]:
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(esave,origin='lower',aspect='auto')
plt.subplot(142)
fftex=sf.fftshift(sf.fft2(esave[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(kaxis,np.sqrt(1+3*kaxis**2*vtc**2),'w--')
plt.axis([kmin/4,kmax/4,wmin/4,wmax/4])
#plt.clim(0,500)
plt.subplot(143)
fftex=sf.fftshift(sf.fft(esave,axis=0))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,xx[nx-1],wmin,wmax]))
plt.axis([0,xx[nx-1],wmin/4,wmax/4])
#plt.clim(0,10)
plt.subplot(144)
fftex=sf.fftshift(sf.fft(esave.T,axis=0),axes=0)
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([0,tt[nt-1],wmin,wmax]))
plt.axis([0,tt[nt-1],wmin/4,wmax/4])

plt.show()