Shape Function

We birefly explain shape functions used in a particle-in-cell (PIC) simulation.

%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
import scipy.fftpack as sf
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5

nx=256
ppc=1
nop=ppc*nx
dx=1
lx=dx*nx
xp=lx*np.random.rand(nop)

Nearest Grid Point (NGP)

The first choice is called the nearest grid point method.
This is the most simple one. It can be defined as follows: $${S_i}({x_p}) = \begin{array}{*{20}{c}} {1,\,\,\,\,if\,\,\left| {\frac{{{x_p} - {X_i}}}{{\Delta x}}} \right| \leqslant \frac{1}{2}} \\ {0,\,\,\,\,else} \end{array}$$

%%time
ds=np.zeros(nx)
round=lambda x:(x*2+1)//2
for ip in range(nop):
    ixm=int(round(xp[ip]/dx))
    wxm=1

    if ixm>nx-1: ixm=ixm-nx

    ds[ixm]=ds[ixm]+wxm/ppc

plt.figure(figsize=(12,10))
plt.subplot(211)
plt.plot(ds,'.-')
plt.xlabel('x'); plt.ylabel('density')
plt.subplot(212)
fdsngp=sf.fftshift(sf.fft(ds-1))
plt.plot(abs(fdsngp))
plt.xlabel('k'); plt.ylabel('power')
plt.show()

CPU times: user 453 ms, sys: 31.2 ms, total: 484 ms
Wall time: 704 ms

Cloud in Cell (CIC)

The second example is the cloud in cell (CIC) method.
This is the most popular one. It can be defined as follows: $${S_i}({x_p}) = \begin{array}{*{20}{c}} {1 - \left| {\frac{{{x_p} - {X_i}}}{{\Delta x}}} \right|,\,\,\,\,if\,\,\left| {\frac{{{x_p} - {X_i}}}{{\Delta x}}} \right| \leqslant 1} \\ {0,\,\,\,\,else} \end{array}$$

%%time
ds=np.zeros(nx)
for ip in range(nop):
    ixm=mt.floor(xp[ip]/dx); ixp=ixm+1
    wxp=xp[ip]/dx-ixm;       wxm=1-wxp

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

    ds[ixm]=ds[ixm]+wxm/ppc
    ds[ixp]=ds[ixp]+wxp/ppc 

plt.figure(figsize=(12,10))
plt.subplot(211)
plt.plot(ds,'.-')
plt.xlabel('x'); plt.ylabel('density')
plt.subplot(212)
fdscic=sf.fftshift(sf.fft(ds-1))
plt.plot(abs(fdscic))
plt.xlabel('k'); plt.ylabel('power')
plt.show()

CPU times: user 422 ms, sys: 31.2 ms, total: 453 ms
Wall time: 456 ms

Second-order Spline

The last one is the second-order spline method.
This is recently used to reduce . It can be defined as follows: \begin{gathered} {S_i}({x_p}) = \begin{array}{*{20}{c}} {\frac{3}{4} - {{\left( {\frac{{{x_p} - {X_i}}}{{\Delta x}}} \right)}^2},{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} if{\mkern 1mu} {\mkern 1mu} \left| {\frac{{{x_p} - {X_i}}}{{\Delta x}}} \right|\leq\frac{1}{2}} \\ {0,{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} else} \end{array} \hfill \\ \\ {S_{i \pm 1}}({x_p}) = \begin{array}{*{20}{c}} {\frac{1}{2}{{\left( {\frac{1}{2} \pm \frac{{{x_p} - {X_i}}}{{\Delta x}}} \right)}^2},{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} if{\mkern 1mu} {\mkern 1mu} \left| {\frac{{{x_p} - {X_i}}}{{\Delta x}}} \right|\leq\frac{1}{2}} \\ {0,{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} else} \end{array} \hfill \\ \hfill \\ \end{gathered}

%%time
ds=np.zeros(nx)
round=lambda x:(x*2+1)//2
for ip in range(nop):
    ixx=int(round(xp[ip]/dx)); wx=(xp[ip]/dx-ixx)

    wxx=3./4.-wx**2
    wxm=0.5*(0.5-wx)**2
    wxp=0.5*(0.5+wx)**2
    if ixx>nx-1: ixx=ixx-nx
    ixp=ixx+1; ixm=ixx-1
    if ixp>nx-1: ixp=ixp-nx
    if ixm<0   : ixm=ixm+nx

    ds[ixm]=ds[ixm]+wxm/ppc
    ds[ixx]=ds[ixx]+wxx/ppc
    ds[ixp]=ds[ixp]+wxp/ppc 

plt.figure(figsize=(12,10))
plt.subplot(211)
plt.plot(ds,'.-')
plt.xlabel('x'); plt.ylabel('density')
plt.subplot(212)
fdsspln=sf.fftshift(sf.fft(ds-1))
plt.plot(abs(fdsspln))
plt.xlabel('k'); plt.ylabel('power')
plt.show()

CPU times: user 453 ms, sys: 15.6 ms, total: 469 ms
Wall time: 467 ms

Comparison

Here is the plot of a Fourier spectrum for the three cases.
We can see that the 2nd spline method can reduce short wave length noises.

plt.figure(figsize=(20,6))
plt.plot(abs(fdsngp),'--k', label="NGP")
plt.plot(abs(fdscic),'-b', label="CIC")
plt.plot(abs(fdsspln),'-r', label="2nd spline")
plt.xlabel('k'); plt.ylabel('power'); plt.legend() 
plt.show()