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

\begin{gathered} {\mathbf{v}}_i^{n + 1/2} = {\mathbf{v}}_i^{n - 1/2} + \frac{{{q_i}}}{{{m_i}}}\left( {{{\mathbf{E}}^n} + \frac{{{\mathbf{v}}_i^n}}{c} \times {{\mathbf{B}}^n}} \right)\Delta t, \hfill \\ {\mathbf{x}}_i^{n + 1} = {\mathbf{x}}_i^n + {\mathbf{v}}_i^{n + 1/2}\Delta t, \hfill \\ {{\mathbf{B}}^{n + 1/2}} = {{\mathbf{B}}^n} - \frac{{c\Delta t}}{2}\nabla \times {{\mathbf{E}}^n}, \hfill \\ {{\mathbf{E}}^{n + 1/2}} = - {\mathbf{V}}_i^{n + 1/2} \times {{\mathbf{B}}^{n + 1/2}} - \frac{{\nabla {P_e}^{n + 1/2}}}{{{q_i}N_i^{n + 1/2}}} - \frac{{{{\mathbf{B}}^{n + 1/2}} \times \left( {\nabla \times {{\mathbf{B}}^{n + 1/2}}} \right)}}{{4\pi {q_i}N_i^{n + 1/2}}} \hfill \\ \end{gathered}

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(lx,up,vp)                    #x_n+1
        ds2,ux2,uy2,uz2=accumulate(nop,nx,dx,up,vp)               #N_n+1,V_n+3/4 
        
        ds=0.5*(ds1+ds2)/ppc  #N_n+1/2
        ux=0.5*(ux1+ux2)/ppc  #Vx_n+1/2
        uy=0.5*(uy1+uy2)/ppc  #Vy_n+1/2
        uz=0.5*(uz1+uz2)/ppc  #Vz_n+1/2
        
        efdm[:,:]=ef          #E_n
        
        bf =Faradayslaw(dth,dx,bf,ef)        #B_n+1/2
        ef =Ohmslaw(Te,dx,bf,ds,ux,uy,uz)    #E_n+1/2
        
        bfdm[:,:]=bf                         #B_n+1/2
        
        ###### predictor ###
        efdm=-efdm+2*ef                      #E'_n+1 
        bfdm =Faradayslaw(dth,dx,bfdm,efdm)  #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(lx,updm,vpdm)                         #x'_n+2
        ds2,ux2,uy2,uz2=accumulate(nop,nx,dx,updm,vpdm)                    #N'_n+2,V'_n+7/4 
        
        ds=0.5*(ds1+ds2)/ppc   #N'_n+3/2
        ux=0.5*(ux1+ux2)/ppc   #Vx'_n+3/2
        uy=0.5*(uy1+uy2)/ppc   #Vy'_n+3/2
        uz=0.5*(uz1+uz2)/ppc   #Vz'_n+3/2
        
        bfdm =Faradayslaw(dth,dx,bfdm,efdm)     #B'_n+3/2
        efdm =Ohmslaw(Te,dx,bfdm,ds,ux,uy,uz)   #E'_n+3/2
        
        ef=0.5*(efdm+ef)                        #E_n+1
        bf =Faradayslaw(dth,dx,bf,ef)           #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

@jit('f8[:](f8,f8[:],f8[:,:])')
def particlepush(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):
    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):
    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)
100%|██████████| 2048/2048 [00:52<00:00, 39.20it/s]
CPU times: user 46 s, sys: 1.98 s, total: 47.9 s
Wall time: 52.3 s

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 [7]:
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()