Vlasov Simulation: CIP Scheme

Reference: "Astrophysical simulation technical manual", Chap. 7, T. Amano (in Japanese)

We show a numerical method to solve the Vlasov-Maxwell system. Equations solved are described as follows: \begin{gathered} \frac{{\partial {f_s}}}{{\partial t}} + {\mathbf{v}} \cdot \frac{{\partial {f_s}}}{{\partial {\mathbf{x}}}} + {{\mathbf{F}}_s} \cdot \frac{{\partial {f_s}}}{{\partial {\mathbf{v}}}} = 0, \hfill \\ {{\mathbf{F}}_s} = \frac{{{q_s}}}{{{m_s}}}\left( {{\mathbf{E}} + \frac{{\mathbf{v}}}{c} \times {\mathbf{B}}} \right), \hfill \\ \frac{1}{c}\frac{{\partial {\mathbf{E}}}}{{\partial t}} = \nabla \times {\mathbf{B}} - \frac{{4\pi }}{c}{\mathbf{J}}, \hfill \\ \frac{1}{c}\frac{{\partial {\mathbf{B}}}}{{\partial t}} = - \nabla \times {\mathbf{E}}, \hfill \\ \nabla \cdot {\mathbf{E}} = 4\pi \rho , \hfill \\ \nabla \cdot {\mathbf{B}} = 0, \hfill \\ \rho = \sum\limits_s {\int {{f_s}({\mathbf{x}},{\mathbf{v}})d{\mathbf{v}}} } , \hfill \\ {\mathbf{J}} = \sum\limits_s {\int {{\mathbf{v}}{f_s}({\mathbf{x}},{\mathbf{v}})d{\mathbf{v}}} } \hfill \\ \end{gathered}

The subscript $s$ indicates the type of particle (ion or electron).

1D1V Vlasov-Poisson Simulation

The background magnetic field is ignored in this case. The govering equations are: \begin{gathered} \frac{{\partial {f_s}}}{{\partial t}} + v\frac{{\partial {f_s}}}{{\partial x}} + \frac{{{q_s}}}{{{m_s}}}E\frac{{\partial {f_s}}}{{\partial v}} = 0, \hfill \\ {\nabla ^2}\phi = - 4\pi \rho, \hfill \\ E = - \frac{{\partial \phi }}{{\partial x}}. \hfill \\ \end{gathered}

We employ the spliting scheme proposed by [Cheng&Knorr1976] here. The above first equation are splitted into two parts and solved alternately. As a simplification, ions are fixed. The splitted equations are \begin{gathered} {\mathbf{J}} = \sum\limits_s {\int {{\mathbf{v}}{f_s}({\mathbf{x}},{\mathbf{v}})d{\mathbf{v}}} }, \hfill \\ \frac{{\partial f}}{{\partial t}} + v\frac{{\partial f}}{{\partial x}} = 0, \hfill \\ \frac{{\partial f}}{{\partial t}} - \frac{e}{m}E\frac{{\partial f}}{{\partial v}} = 0. \hfill \\ \end{gathered}

Cheng&Knorr(1976) shows that when you solve these as follows, the accuracy of this scheme becomes second order, \begin{gathered} {f^*}\left( {x,v} \right) = {f^n}\left( {x - v\Delta t/2,v} \right), \hfill \\ {f^{**}}\left( {x,v} \right) = {f^*}\left( {x,v + \frac{e}{m}{E^*}\left( x \right)\Delta t} \right), \hfill \\ {f^{n + 1}}\left( {x,v} \right) = {f^{**}}\left( {x - v\Delta t/2,v} \right). \hfill \\ \end{gathered} .

CIP Method

We need to choose a method to advance the distribution function. The CIP method is a high order accurate solver for the linear advection equation; $$\frac{{\partial f}}{{\partial t}} + c\frac{{\partial f}}{{\partial x}} = 0$$.

In this method, values between $(i,i+1)$ are interpolated with a cubic function, $$f\left( x \right) = a{X^3} + b{X^2} + cX + d,\,\,X = x - {x_i}$$. The coefficients $a$, $b$, $c$, and $d$ are determined by \begin{gathered} a = \frac{{{g_i} + {g_{i + 1}}}}{{\Delta {x^2}}} + \frac{{2\left( {{f_i^n} - {f_{i + 1}^n}} \right)}}{{\Delta {x^3}}} \hfill \\ b = \frac{{3\left( {{f_{i + 1}^n} - {f_i^n}} \right)}}{{\Delta {x^2}}} - \frac{{2{g_i} + {g_{i + 1}}}}{{\Delta x}} \hfill \\ c = {g_i} \hfill \\ d = {f_i}. \hfill \\ \end{gathered} Here, $g_i$ is a drivative in space at $i$.

After $\Delta t$, $f$ and $g$ become \begin{gathered} {f^{n + 1}_i} = a{\xi ^3} + b{\xi ^2} + c\xi + d \hfill \\ {g^{n + 1}_i} = 3a{\xi ^2} + 2b\xi + c \hfill \\ \end{gathered} where $\xi$ is $c\Delta t$.

Advection equation: CIP method

Here, we show an example of a simulation solved by CIP method. We solve the advection equation. $$\frac{{\partial f}}{{\partial t}} + c\frac{{\partial f}}{{\partial x}} = 0$$.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
def CIP(dx,nx,dt,nt):
    c=1
    #dx=-dx
    f =np.r_[np.zeros(int(nx/4)), np.ones(int(nx/2)), np.zeros(int(nx/4))]
    df=(np.roll(f,1)-np.roll(f,-1))/(2*dx)
    xi=c*dt
    fhst=np.zeros((nt,nx))
    fhst[0,:]=f    
    for it in range(1,nt):
        a=(df+np.roll(df,1))/dx**2+2*(f-np.roll(f,1))/dx**3
        b=3*(np.roll(f,1)-f)/dx**2-(2*df+np.roll(df,1))/dx
        c=df
        d=f
        
        f =a*xi**3+b*xi**2+c*xi+d
        df=3*a*xi**2+2*b*xi+c
        
        fhst[it,:]=f
    
    plt.figure(figsize=(12,4))
    plt.subplot(121);plt.plot(fhst[0,:],'k--'); plt.plot(fhst[nt-1,:])
    plt.title('initial and last state')
    plt.subplot(122);plt.pcolor(fhst)
    plt.show()
In [3]:
%%time
CIP(1,128,0.1,1280)

Wall time: 6.1 s

Application of CIP Method to Vlasov-Poison System

We need to find $\partial f^*/\partial v$.
Nakamura&Yabe(1999) suggests to use \begin{gathered} \frac{{\partial f_{ix,iv}^*}}{{\partial v}} = \frac{{\partial {f_{ix,iv}}}}{{\partial v}} - \frac{{\Delta t}}{{4\Delta v}}\left[ {{v_{iv + 1}}\left( {\frac{{\partial f_{ix,iv + 1}^*}}{{\partial x}} + \frac{{\partial {f_{ix,iv + 1}}}}{{\partial x}}} \right) - {v_{iv - 1}}\left( {\frac{{\partial f_{ix,iv - 1}^*}}{{\partial x}} + \frac{{\partial {f_{ix,iv - 1}}}}{{\partial x}}} \right)} \right] \hfill \\ \frac{{\partial f_{ix,iv}^{**}}}{{\partial x}} = \frac{{\partial f_{ix,iv}^*}}{{\partial v}} + \frac{{\Delta t}}{{4\Delta x}}\left[ {\frac{e}{m}{E_{ix + 1}}\left( {\frac{{\partial f_{ix + 1,iv}^{**}}}{{\partial v}} + \frac{{\partial f_{ix + 1,iv}^*}}{{\partial v}}} \right) - \frac{e}{m}{E_{i - 1}}\left( {\frac{{\partial f_{ix - 1,iv}^{**}}}{{\partial v}} + \frac{{\partial f_{ix - 1,iv}^*}}{{\partial v}}} \right)} \right] \hfill \\ \end{gathered}

Summary of the Simulation

We write a time advance by CIP method with $\Delta t$ as \begin{gathered} \left( {{f^*},\frac{{\partial {f^*}}}{{\partial x}}} \right) &=& {\text{CIP}}\left( {f,\frac{{\partial f}}{{\partial x}},v,\Delta t} \right) \hfill \\ \frac{{\partial {f^*}}}{{\partial v}} &=& {\text{FDM}}\left( {\frac{{\partial f}}{{\partial v}},\frac{{\partial f}}{{\partial x}},\frac{{\partial {f^*}}}{{\partial x}},v,\Delta t} \right). \hfill \\ \end{gathered}

Then, the splitting method can be written as

  1. step1 \begin{gathered} \left( {{f^*},\frac{{\partial {f^*}}}{{\partial x}}} \right) &=& {\text{CIP}}\left( {{f^n},\frac{{\partial {f^n}}}{{\partial x}},v,\frac{{\Delta t}}{2}} \right) \hfill \\ \frac{{\partial {f^*}}}{{\partial v}} &=& {\text{FDM}}\left( {\frac{{\partial {f^n}}}{{\partial v}},\frac{{\partial {f^n}}}{{\partial x}},\frac{{\partial {f^*}}}{{\partial x}},v,\frac{{\Delta t}}{2}} \right) \hfill \\ \end{gathered}
  2. step2 \begin{gathered} \left( {{f^{**}},\frac{{\partial {f^{**}}}}{{\partial v}}} \right) &=& {\text{CIP}}\left( {{f^*},\frac{{\partial {f^*}}}{{\partial v}},{E^*},\Delta t} \right) \hfill \\ \frac{{\partial {f^{**}}}}{{\partial x}} &=& {\text{FDM}}\left( {\frac{{\partial {f^*}}}{{\partial x}},\frac{{\partial {f^*}}}{{\partial v}},\frac{{\partial {f^{**}}}}{{\partial x}},{E^*},\Delta t} \right) \hfill \\ \end{gathered}
  3. step3 \begin{gathered} \left( {{f^{n + 1}},\frac{{\partial {f^{n + 1}}}}{{\partial x}}} \right) &=& {\text{CIP}}\left( {{f^{**}},\frac{{\partial {f^{**}}}}{{\partial x}},v,\frac{{\Delta t}}{2}} \right) \hfill \\ \frac{{\partial {f^{n + 1}}}}{{\partial v}} &=& {\text{FDM}}\left( {\frac{{\partial {f^{**}}}}{{\partial v}},\frac{{\partial {f^{**}}}}{{\partial x}},\frac{{\partial {f^{n + 1}}}}{{\partial x}},v,\frac{{\Delta t}}{2}} \right) \hfill \\ \end{gathered}

Alongside these steps, we have to solve the Poisson equation and obtain $E^*_x$.

Test Problem 1: Linear Landau Damping

You can see this test problem in these articles ([Cheng&Knorr1976], [Nakamura&Yabe1999], and [Filbert+2001]).

The initical distribution function is $${f_0}\left( {x,v} \right) = \frac{1}{{\sqrt 2 }}\exp \left( { - \frac{{{v^2}}}{2}} \right)\left( {1 + A\cos kx} \right),$$ where $A=0.01$ and $k=0.5$. The system size for real space is $4\pi$ and $-4<v<4$ for velocity space.

From the linear analysis of the problem, a damping rate is $\gamma\sim0.1533$. One can see that the result of the Vlasov simulation corresponds to the theory.

In [22]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
import matplotlib.animation as animation
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [18]:
nx=32; nt=500
dt=1/8; gam=0.1533
tt=np.arange(nt)*dt
f    =open('exhst.d'  ,'rb')
exhst=np.fromfile(f, dtype='float64').reshape(-1,nx)
plt.semilogy(tt,abs((fftexhst[:,15])))
Ek=0.3*np.exp(-gam*tt)
plt.plot(tt,Ek,'k--')
plt.xlabel(r'$t\omega_{pe}$'); plt.ylabel(r'$|E_k(k=0.5)|$')
plt.show()

Test Problem 2: Two Stream Instability

The initial distribution function is $${f_0}\left( {x,v} \right) = \frac{1}{{\sqrt 2 }}{v^2}\exp \left( { - \frac{{{v^2}}}{2}} \right)\left( {1 - A\cos kx} \right),$$ where A=0.05 and k=0.5. The system size is $4\pi$ and the cutoff velocity is $5$.

Animation

In [29]:
nx=32; nt=500
dt=1/8
plt.rcParams['animation.html'] = 'html5'
#fig,(ax1,ax2)=plt.subplots(1,3,figsize=(10,5))
fig=plt.figure(figsize=(15,5))
ax = fig.add_subplot(111)

def update_anim(i):
    it =(i+1)*5
    f  =open('data/f%05d.d' %(it),'rb') 
    fnc=np.fromfile(f,dtype='float64').reshape(-1,nx)
    
    ax.clear()

    ax.imshow(fnc,origin='lower',aspect='auto',cmap='gist_heat',extent=[0,4*np.pi,-4,4])
    #ax.imshow(rho,origin='lower',aspect='auto',cmap='gist_heat')
    ax.set_title(r"$t\omega_{pe}=%.2f$" %(it*dt))
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$v$')
    
anim=animation.FuncAnimation(fig,update_anim,frames=100)    
plt.close()
anim
Out[29]:

Time Evolution of the energy of the electric field

In [20]:
tt=np.arange(nt)*dt
f=open('data/exhst.d','rb')
exhst=np.fromfile(f, dtype='float64').reshape(-1,nx)
#plt.imshow(exhst,aspect='auto')
plt.plot(tt,np.mean(exhst**2,axis=1))
plt.xlabel(r'$t\omega_{pe}$'); plt.ylabel(r'$Electric Field Energy$')
plt.show()