Hybrid Simulation¶
Reference: D. Winske, 2003, Space Plasma Simulation. Edited by J. Büchner, C. Dum, and M. Scholer., Lecture Notes in Physics, vol. 615, p.136-165
Basic Assumptions and Equations¶
\begin{gathered} {m_e} \to 0 \hfill \\ e{n_e} = {q_i}{m_i} \hfill \\ \frac{{\partial {\mathbf{E}}}}{{\partial t}} \to 0 \hfill \\ {m_i}\frac{{d{{\mathbf{v}}_i}}}{{dt}} = {q_i}\left( {{\mathbf{E}} + \frac{{{{\mathbf{v}}_i}}}{c} \times {\mathbf{B}}} \right), \hfill \\ \frac{{d{{\mathbf{x}}_i}}}{{dt}} = {{\mathbf{v}}_i}, \hfill \\ {\mathbf{E}} = - {{\mathbf{V}}_i} \times {\mathbf{B}} - \frac{{\nabla {P_e}}}{{{q_i}{N_i}}} - \frac{{{\mathbf{B}} \times \left( {\nabla \times {\mathbf{B}}} \right)}}{{4\pi {q_i}{N_i}}}, \hfill \\ \nabla \times {\mathbf{B}} = \frac{{4\pi }}{c}{\mathbf{J}}, \hfill \\ \frac{{\partial {\mathbf{B}}}}{{\partial t}} = - c\nabla \times {\mathbf{E}}. \hfill \\ \end{gathered}Discretization¶
- Particle positions and velocites are leap-frogged in time
- Electromagnetic fields, density, and current are defined on co-locatted grids
Perdictor-Corrector¶
We need a scheme to obtain ${\bf E}_{n+1}$.
Step 1: Predicting ${\bf E}^\prime_{n+1}$ and ${\bf B}^\prime_{n+1}$ \begin{gathered} {{\mathbf{E}}^{\prime n + 1}} = - {{\mathbf{E}}^n} + 2{{\mathbf{E}}^{n + 1/2}}, \hfill \\ {{\mathbf{B}}^{\prime n + 1}} = {{\mathbf{B}}^{n + 1/2}} - \frac{{c\Delta t}}{2}\nabla \times {{\mathbf{E}}^{\prime n + 1}} \hfill \\ \end{gathered}
Step 2: Advancing particles through the predicted fields to obtain ${\bf V}_i^{\prime n+3/2}$ and $N_i^{\prime n+3/2}$
Step 3: Advancing the predicted fields to $n+3/2$ $${{\mathbf{B}}^{\prime n + 3/2}} = {{\mathbf{B}}^{\prime n + 1}} - \frac{{c\Delta t}}{2}\nabla \times {{\mathbf{E}}^{\prime n + 1}}$$ ${\bf E}^{\prime n+3/2}$ is otained from the Ohm's law
Step 4: Correcting fields \begin{gathered} {{\mathbf{E}}^{n + 1}} = \frac{1}{2}\left( {{{\mathbf{E}}^{n + 1/2}} + {{\mathbf{E}}^{\prime n + 3/2}}} \right), \hfill \\ {{\mathbf{B}}^{n + 1}} = {{\mathbf{B}}^{n + 1/2}} - \frac{{c\Delta t}}{2}\nabla \times {{\mathbf{E}}^{n + 1}} \hfill \\ \end{gathered}
Other Schemes¶
- The Bunemann-Boriss method is applied for a particle acceleratio
- Faraday's law is advanced by FTCS (Forward Time Center Scheme)
Code¶
In [1]:
%config InlineBackend.figure_format = 'retina'
from numba import jit,f8,u8
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,fftshift
from tqdm import tqdm
from IPython import display
In [2]:
def hybr1d(vth,ppc,dx,nx,dt,nt):
#global jp,jm
nop=ppc*nx
Te=0.1
lx=dx*nx; dth=dt/2
q=1; m=1
qomdt=q*dt/(2*m)
xx=dx*np.arange(nx)
jj=np.arange(nx); jp=np.r_[np.arange(nx-1)+1, 0]; jm=np.r_[nx-1, np.arange(nx-1)]
up=np.linspace(0,lx-lx/nop,nop)
vp=vth*np.random.randn(3,nop)
bf =np.zeros((3,nx)); ef =np.zeros((3,nx))
b0x=1; b0z=0.
bf[0,:]=b0x; bf[2,:]=b0z
bfdm =np.zeros((3,nx))
efdm =np.zeros((3,nx))
updm =np.zeros(nop)
vpdm =np.zeros((3,nop))
byhst=np.zeros((nt,nx))
enhst=np.zeros(nt)
for it in tqdm(range(nt)):
vp =particleacc(nop,nx,dx,qomdt,up,vp,bf,ef) #v_n+1/2
ds1,ux1,uy1,uz1=accumulate(nop,nx,dx,up,vp) #N_n, V_n+1/4
up =particlepush(dt,lx,up,vp) #x_n+1
ds2,ux2,uy2,uz2=accumulate(nop,nx,dx,up,vp) #N_n+1,V_n+3/4
ds1=smooth(ds1,jp,jm); ds2=smooth(ds2,jp,jm)
ux1=smooth(ux1,jp,jm); ux2=smooth(ux2,jp,jm)
uy1=smooth(uy1,jp,jm); uy2=smooth(uy2,jp,jm)
uz1=smooth(uz1,jp,jm); uz2=smooth(uz2,jp,jm)
ds=0.5*(ds1+ds2)/ppc #N_n+1/2
ux=0.5*(ux1+ux2)/ppc/ds #Vx_n+1/2
uy=0.5*(uy1+uy2)/ppc/ds #Vy_n+1/2
uz=0.5*(uz1+uz2)/ppc/ds #Vz_n+1/2
efdm[:,:]=np.copy(ef) #E_n
bf =Faradayslaw(dth,dx,bf,ef,jp,jm) #B_n+1/2
ef =Ohmslaw(Te,dx,bf,ds,ux,uy,uz,jp,jm) #E_n+1/2
bfdm[:,:]=np.copy(bf) #B_n+1/2
###### predictor ###
efdm=-efdm+2*ef #E'_n+1
bfdm =Faradayslaw(dth,dx,bfdm,efdm,jp,jm) #B'_n+1
###### corrector ###
updm[:]=up; vpdm[:,:]=vp #x'_n+1, v'_n+1/2
vpdm =particleacc(nop,nx,dx,qomdt,updm,vpdm,bfdm,efdm) #V'_n+1/2
ds1,ux1,uy1,uz1=accumulate(nop,nx,dx,updm,vpdm) #N'_n+1, V'_n+5/4
updm =particlepush(dt,lx,updm,vpdm) #x'_n+2
ds2,ux2,uy2,uz2=accumulate(nop,nx,dx,updm,vpdm) #N'_n+2,V'_n+7/4
ds1=smooth(ds1,jp,jm); ds2=smooth(ds2,jp,jm)
ux1=smooth(ux1,jp,jm); ux2=smooth(ux2,jp,jm)
uy1=smooth(uy1,jp,jm); uy2=smooth(uy2,jp,jm)
uz1=smooth(uz1,jp,jm); uz2=smooth(uz2,jp,jm)
ds=0.5*(ds1+ds2)/ppc #N'_n+3/2
ux=0.5*(ux1+ux2)/ppc/ds #Vx'_n+3/2
uy=0.5*(uy1+uy2)/ppc/ds #Vy'_n+3/2
uz=0.5*(uz1+uz2)/ppc/ds #Vz'_n+3/2
bfdm =Faradayslaw(dth,dx,bfdm,efdm,jp,jm) #B'_n+3/2
efdm =Ohmslaw(Te,dx,bfdm,ds,ux,uy,uz,jp,jm) #E'_n+3/2
ef=0.5*(efdm+ef) #E_n+1
bf =Faradayslaw(dth,dx,bf,ef,jp,jm) #B_n+1
byhst[it,:]=bf[2,:]
enhst[it] =np.mean(0.5*(vp[0,:]**2+vp[1,:]**2+vp[2,:]**2))
#display.clear_output(True)
#plt.figure(figsize=(18,6))
#plt.plot(up,vp[0,:],'.')
#plt.show()
return locals()
@jit('f8[:,:](u8,u8,f8,f8,f8[:],f8[:,:],f8[:,:],f8[:,:])',nopython=True)
def particleacc(nop,nx,dx,qomdt,up,vp,bf,ef):
for ip in range(nop):
ixm=int(up[ip]/dx)
ixp=ixm+1
wxp=up[ip]/dx-ixm
wxm=1-wxp
if ixp>nx-1: ixp=ixp-nx
bbx=qomdt*(wxm*bf[0,ixm]+wxp*bf[0,ixp])
bby=qomdt*(wxm*bf[1,ixm]+wxp*bf[1,ixp])
bbz=qomdt*(wxm*bf[2,ixm]+wxp*bf[2,ixp])
eex=qomdt*(wxm*ef[0,ixm]+wxp*ef[0,ixp])
eey=qomdt*(wxm*ef[1,ixm]+wxp*ef[1,ixp])
eez=qomdt*(wxm*ef[2,ixm]+wxp*ef[2,ixp])
vmx=vp[0,ip]+eex; vmy=vp[1,ip]+eey; vmz=vp[2,ip]+eez
v0x=vmx+(vmy*bbz-vmz*bby)
v0y=vmy+(vmz*bbx-vmx*bbz)
v0z=vmz+(vmx*bby-vmy*bbx)
tt =2/(1+bbx**2+bby**2+bbz**2)
vpx=vmx+tt*(v0y*bbz-v0z*bby)
vpy=vmy+tt*(v0z*bbx-v0x*bbz)
vpz=vmz+tt*(v0x*bby-v0y*bbx)
vp[0,ip]=vpx+eex; vp[1,ip]=vpy+eey; vp[2,ip]=vpz+eez
return vp
@jit('(u8,u8,f8,f8[:],f8[:,:])')
def accumulate(nop,nx,dx,up,vp):
ds=np.zeros(nx)
ux=np.zeros(nx); uy=np.zeros(nx); uz=np.zeros(nx)
for ip in range(nop):
ixm=int(up[ip]/dx)
ixp=ixm+1
wxp=up[ip]/dx-ixm
wxm=1-wxp
if ixp>nx-1: ixp=ixp-nx
ds[ixm]=ds[ixm]+wxm
ds[ixp]=ds[ixp]+wxp
ux[ixm]=ux[ixm]+wxm*vp[0,ip]
ux[ixp]=ux[ixp]+wxp*vp[0,ip]
uy[ixm]=uy[ixm]+wxm*vp[1,ip]
uy[ixp]=uy[ixp]+wxp*vp[1,ip]
uz[ixm]=uz[ixm]+wxm*vp[2,ip]
uz[ixp]=uz[ixp]+wxp*vp[2,ip]
#### smopthing ###
#ds=ds[jp]/4+ds/2+ds[jm]/4
#ux=ux[jp]/4+ux/2+ux[jm]/4
#uy=uy[jp]/4+uy/2+uy[jm]/4
#uz=uz[jp]/4+uz/2+uz[jm]/4
return ds,ux,uy,uz
def smooth(zz,jp,jm):
zz=zz[jp]/4+zz/2+zz[jm]/4
return zz
@jit('f8[:](f8,f8,f8[:],f8[:,:])')
def particlepush(dt,lx,up,vp):
up=up+dt*vp[0,:]
up[up>lx]=up[up>lx]-lx
up[up<0 ]=up[up<0 ]+lx
return up
#@jit('f8[:,:](f8,f8,f8[:,:],f8[:,:])')
def Faradayslaw(dth,dx,bf,ef,jp,jm):
bf[1,:]=bf[1,:]+0.5*dth/dx*(ef[2,jp]-ef[2,jm])
bf[2,:]=bf[2,:]-0.5*dth/dx*(ef[1,jp]-ef[1,jm])
return bf
#@jit('f8[:,:](f8,f8,f8[:,:],f8[:],f8[:],f8[:],f8[:])')
def Ohmslaw(Te,dx,bf,ds,ux,uy,uz,jp,jm):
ef=np.zeros((3,nx))
bp2=bf[1,:]**2+bf[2,:]**2
ef[0,:]=-uy*bf[2,:]+uz*bf[1,:]-Te*(ds[jp]-ds[jm])/(2*dx*ds)-(bp2[jp]-bp2[jm])/(4*dx*ds)
ef[1,:]=-uz*bf[0,:]+ux*bf[2,:] +bf[0,:]*(bf[1,jp]-bf[1,jm])/(2*dx*ds)
ef[2,:]=-ux*bf[1,:]+uy*bf[0,:] +bf[0,:]*(bf[2,jp]-bf[2,jm])/(2*dx*ds)
return ef
In [3]:
%%time
vth=0.2; ppc=100; dx=0.5; nx=1024; dt=0.1; nt=2048
data=hybr1d(vth,ppc,dx,nx,dt,nt)
locals().update(data)
Simulation Results: Linear Waves¶
Time Evolution of Particle Energy¶
In [4]:
tt=np.arange(nt)*dt
plt.plot(tt,enhst)
plt.xlabel(r'$t\Omega_p$')
plt.ylabel(r'Averaged Energy')
plt.show()
Time-Space Evolution of $B_y$¶
In [5]:
plt.imshow(byhst,origin='lower',aspect='auto',cmap='jet',extent=[0,nx*dx,0,nt*dt])
plt.colorbar()
plt.xlabel(r'$x/(v_A/\Omega_p)$')
plt.ylabel(r'$t\Omega_p$')
plt.title(r'$B_y$')
plt.show()
Dispersion Relation¶
In [6]:
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)
fftby=fftshift(fft2(byhst[:,::-1]))
plt.imshow(abs(fftby),origin='lower',aspect='auto',extent=[kmin,kmax,wmin,wmax])
plt.xlabel(r'$k(v_A/\Omega_p)$')
plt.ylabel(r'$\omega/\Omega_p$')
plt.colorbar()
plt.clim(0,100)
plt.xlim(kmin/4, kmax/4)
plt.ylim(wmin/16,wmax/16)
plt.show()
