Chapt 7: Particle-In-Cell

Author
Affiliation

Feng Wan

B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaoton University

Show the code
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from IPython.display import Markdown
from sympy.solvers.ode.systems import dsolve_system
from sympy import *
plt.rcParams.update({
    "text.usetex": True,
    "font.family": "Helvetica"
})

1 激光等离子体相互作用

1.1 激光等离子体相互作用

Laser plasma phenomenon and applications:

  • \(I \lesssim 10^{15} \mathrm{W/cm^2}\) 多光子电离
  • \(I \gtrsim 3.51\times 10^{16} \mathrm{W/cm^2}\) 完全电离
  • \(I \sim ~10^{19} \mathrm{W/cm^2}\) 激光等离子体加速, ICF等
  • \(I \sim 10^{21} \mathrm{W/cm^2}\) 辐射压力加速等

1.2 激光等离子体相互作用

Quantum parameters controlling the nonlinear Compton scattering and pair conversion: \[\eta (\chi_e) = \frac{e\hbar}{m^3c^4}|F_{\mu \nu}p^\nu| \sim \gamma E_\perp / E_\mathrm{cr}\] \[ \chi (\chi_\gamma) = \frac{e\hbar}{ m^3c^4}|F_{\mu \nu}k^\nu| \sim \frac{\hbar \omega}{m_e c^2} E_\perp / E_\mathrm{cr} \]

(The ratio between field of particle experienced and the Schwinger field)


Assumption:

  • time variation \(\gg\) coherent interaction time: \(t_\mathrm{coh} \simeq (E_\mathrm{cr}/E_0) / (\hbar/mc^2) \Longrightarrow a \gg 1\): uniform static field

  • weak field approximation: \[ \begin{eqnarray}f &=& |E^2 - B^2|/E_\mathrm{cr}^2 \ll 1, \\ g &=& |{\bf E} \cdot {\bf B}|\frac{e\hbar}{m^3c^4}|F_{\mu \nu}p^\nu| \sim \gamma E_\perp / E_\mathrm{cr}; \end{eqnarray}\] \[\eta \gg f, g\] Transition probability only depends on \(\eta, \chi\).

2 粒子模拟方法(particle-in-cell, PIC)

2.1 数值模拟方法

  • 流体描述是从宏观角度出发, 研究等离子体大范围, 长时间的演化性质, 将微观研究得到的吸收系数或输运系数等作为已知的条件, 数值求解磁流体方程。

  • 动理学描述是从微观角度出发, 建立在等离子体粒子组态与速度空间分布以及由这些粒子产生的微观场基础之上的描述。等离子体动理学描述有两种不同的方法:

    • 直接求解Vlasov方程(考虑碰撞效应时为Fokker-Planck方程)

      动理学方程(Vlasov方程或Fokker-Planck方程)的求解由于存在一个多维相空间的分布函数,数值求解时要进行离散化处理,这样容易产生非物理的多束流失真,掩盖真正的物理解。
    • 追踪在自洽场和外加场作用下的带电粒子运动的方法, 称为粒子模拟方法(particle-in-cell, PIC)


无论流体方程还是动理学方程,在建立的时候都对作为统计系统的等离子体作了光滑的近似,抹去了他们固有的统计起伏效应,而这些效应在一定的条件下可以发展成为像湍流这样的重要物理现象。

2.2 数值模拟方法

数值模拟方法

数值模拟方法

2.3 粒子模拟方法

  • 等离子体粒子模拟方法(particle-in-cell)是通过跟踪大量的带电粒子(超粒子,巨粒子)在它们的自洽场和外加电磁场中的运动来研究等离子体集体性质的动理学模拟方法。

  • 从原则上讲,这种方法考虑等离子体运动最齐全,最能反映实际等离子体的运动,在一定意义上,可以代替实验的功能。又因为它是模拟计算,每一个模拟的带电粒子的全部运动都存在计算机里,因而可以提供任何详细的等离子体运动的信息,在这一点上,它又高于实验。现在,等离子体粒子模拟方法已经成为研究等离子体非线性效应的强有力的工具。

  • 等离子体粒子模拟的基本思想是给定初始时刻粒子的位置速度,用平均近似的方法求得空间各处的电荷密度电流密度分布,求解Maxwell方程组得到空间各处的电场和磁场,再用平均近似的方法求出粒子感受到的电场强度和磁场强度,每个粒子所受的洛伦兹力便可求得,从而下一时刻各个粒子的位置和速度就可以通过求解运动方程求出。如此循环求解,就可以知道电磁场和每个粒子运动的演化情况。所有关于场和粒子的基本量都可以用计算机储存起来,任何宏观量和统计量都可以通过这些基本量求出。

2.4 粒子模拟方法

2.5 粒子运动

\[\begin{eqnarray} \frac{\mathbf{p}^{n+1/2} - \mathbf{p}^{n-1/2}}{\Delta t} &=& \frac{q}{m}\left(\mathbf{E}^n + \frac{\mathbf{p}^n}{\gamma^n} \times \mathbf{B}^n\right), \label{lorentz1} \\ \frac{\mathbf{x}^{n+1} - \mathbf{x}^n}{\Delta t} &=& \mathbf{v}^{n+1/2}, \label{lorentz2} \end{eqnarray}\]
\[\begin{equation} \begin{split} \left(\frac{{\rm d}{\bf S}}{{\rm d}t}\right)_{T}= ~&\mathbf{S} \times {\bf \Omega} \equiv {\bf S}\times \frac{1}{\gamma}\left[-\left(\frac{g}{2}-1\right)\frac{\gamma_e}{\gamma_e+1}\left({\bf \beta}\cdot{\bf B}\right)\cdot{\bf \beta} \right.\\ &\left. +\left(\frac{g}{2}-1+\frac{1}{\gamma_e}\right){\bf B}-\left(\frac{g}{2}-\frac{\gamma_e}{\gamma_e+1}\right)\bf{\beta}\times{\bf E}\right], \end{split}\label{bmteqn} \end{equation}\]

2.6 电磁场演化

完整Maxwell方程组 :

\[\begin{eqnarray*} \nabla \cdot \mathbf{D} &=& \rho \\ \nabla \cdot \mathbf{B} &=& 0 \\ \nabla \times \mathbf{E} &=& -\frac{\partial \mathbf{B}}{\partial t} \\ \nabla \times \mathbf{H} &=& \frac{\partial \mathbf{D}}{\partial t} + \mathbf{J} \\ \mathbf{D} &=& \mathbf{E} + \mathbf{P}, \mathbf{B} = \mathbf{H} + \mathbf{M} \end{eqnarray*}\]

完全电离的等离子体:\(\mathbf{P} = \mathbf{M} = 0\)。一般等离子体:\(\mathbf{P} \ll E, \mathbf{M} \ll B\)

\[\begin{eqnarray*} \nabla \cdot \mathbf{E} &=& \rho \\ \nabla \cdot \mathbf{B} &=& 0 \\ \nabla \times \mathbf{E} &=& -\frac{\partial \mathbf{B}}{\partial t} \\ \nabla \times \mathbf{B} &=& \frac{\partial \mathbf{E}}{\partial t} + \mathbf{J}. \end{eqnarray*}\]

2.7 有限粒子大小

有限粒子大小

真实等离子体系统中带电粒子数目巨大,由于计算机容量及速度的限制,直接对单个粒子运算是不实际的。实际模拟时引入超粒子(superparticle)或巨粒子(macroparticle)的概念,即用一个模拟粒子代替分布在空间小范围内的粒子团。

在等离子体空间中一点的周围,每个带点粒子对电磁场的贡献及电磁场对粒子的作用都基本相同,故周围这些大量带电粒子的运动规律基本相同,因而可以不必跟踪计算每一个粒子,只计算代表这些粒子的一个粒子就可以了,即超粒子。这与我们常常对大于德拜长度的无碰撞等离子体的集体效应感兴趣是一致的。

2.8 有限粒子大小

数值模拟方法

数值模拟方法
有限大小粒子

点状超粒子不行(或不理想),因为它将近程碰撞效应放大了,引起了非物理效应(当两团粒子重叠时,它们的库仑作用力会减小,完全重叠时库仑作用力为0. 但当把两团电子气体集中在两个点上时,它们距离远时作用力与团粒子一样,但当点粒子靠近时,库仑力就会趋于无穷)。为解决这一困难,C.K.Birdsall等人提出了有限大小粒子(FSP:finite-size particle)模型。有限大小粒子模型的基本思想是: 模拟粒子是电荷在空间有一定分布、有一定形状的粒子云,粒子在距离近时可以互相穿透甚至重合。

2.9 有限粒子大小

Show the code
x = np.linspace(-2, 2, 100)
def y1(x):
    a = 1
    res = np.zeros(x.shape)
    res[np.abs(x) < a] = 1.0 / a
    return res

def y2(x):
    a = 1
    res = np.zeros(x.shape)
    res = np.exp(-x**2.0 / (2.0 * a**2.0)) / (2.0 * np.sqrt(np.pi) * a)
    return res

str1 = (r"\begin{eqnarray*}"
       r"S(r) &=& \frac{1}{V_n a^n}, ~~ r < a, \\ "
       r"S(r) &=& 0, ~~ r > a \end{eqnarray*}")

str2 = (r"$$S(r) = \frac{\exp\left\{\frac{r^2}{2a_0^2}\right\}}{2\pi^{n/2}a_0^n}$$")

plt.figure(figsize = (10, 3))
plt.subplot(121)
plt.plot(x, y1(x))
plt.text(-1, 1.3, str1, fontsize = 15);

plt.subplot(122)
plt.plot(x, y2(x))
plt.text(-1, 0.32, str2, fontsize = 15);

有限大小粒子

和点粒子不同,有限大小粒子的电荷不再是集中于一点(\(\delta\)函数),而是按一定的形状(称为有限大小粒子形状因子,用\(S(r)\)表示) 分布在它的有限空间里,从而形成了空间电荷密度分布。至于采用何种电荷分布的有限大小粒子,则不受任何限制(取决于weighting方法)。通常所用的有限大小粒子形状因子有球形分布(1D: 矩形分布)、三角形分布、高斯分布等。

2.10 有限粒子大小

2.11 有限粒子大小

  • 形状因子需要满足 \[ \int S(\mathbf{r})d^n \mathbf{r} = 1\]
  • 有限大小粒子形成的空间电荷、电流密度分布 \[ \rho(\mathbf{r}) = \sum_i q_i S(\mathbf{r} - \mathbf{r}_i) \] \[\mathbf{J}(\mathbf{r}) = \sum_i q_i \mathbf{v}_i S(\mathbf{r} - \mathbf{r}_i)\]

2.12 直角坐标系插值

  • 零阶权重法 (zero-order weighting, nearest-grid-point (NGP) )
  • L-E: 计算位于 jth网格点附近±Δx/2范围内的粒子个数,并认为这就是位于jth网格的粒子个数,从而网格点的密度为

\[ n_j = N(j) / \Delta x\]

Zero-th order

Zero-th order

2.13 直角坐标系插值

NGP特点:

  • 密度分布(形状因子)为方形。这种有限大小粒子(particle-in-cell)的引入可以减少近距离碰撞,而不改变等离子体长程相互作用。
  • weighting时只需考虑1个网格点,计算速度快。
  • 方形的密度在边界处不连续,出现初始非电中性的电荷分离,而等离子体中存在强的恢复电中性的倾向,于是引发噪声。这种噪声在很多等离子体过程中不允许。(很少用)

更高阶插值函数

2.14 直角坐标系插值

  • 一阶权重法(线性插值) (first-order weighting, cloud in cell (CIC) 模型 )
  • L-E: 将有限大小粒子看成密度均匀的粒子云,设粒子云的宽度为\(\Delta x\), 中心为\(x_i\) , weighting时将此粒子云上的量分到包含粒子云中心所在的两个网格点上,每个网格点上用NGP方法weighting.

2.15 直角坐标系插值

CIC特点:

  • 密度分布(形状因子)为三角形。 三角形的密度在边界处连续,降低了叠加到等离子体密度和场上的噪声扰动,因而在等离子体模拟中得到了广泛应用。
  • CIC在每个粒子上的计算时间比NGP多,但当给定同一噪声水平时,CIC可以比NGP允许更大的网格和较少的粒子,从而能有效补偿每个粒子上计算时多用的时间。

2.16 直角坐标系插值

  • 零阶:密度分布(形状因子)为方形,不连续。计算速度快。
  • 一阶:密度分布(形状因子)为三角形,连续,但不光滑。比零阶耗时。
  • 高阶:密度分布(形状因子)连续,光滑。(C.K. Birdsall, A. B. Langdon, Plasma Physics Via Computer Simulation, chapter 8-8) 能很好降低噪声,减少非物理效应,但更耗时。

2.17 直角坐标系插值

面积权重法:

\[\begin{eqnarray*} &&q_{j, k} = q_i \frac{A_d}{\Delta x \Delta y} = q_i \frac{(\Delta x - x)(\Delta y - y)}{\Delta x \Delta y} \\ &&q_{j+1, k} = q_i \frac{A_c}{\Delta x \Delta y} = q_i \frac{x(\Delta y - y)}{\Delta x \Delta y} \\ &&q_{j+1, k+1} = q_i \frac{A_b}{\Delta x \Delta y} = q_i \frac{xy}{\Delta x \Delta y} \\ &&q_{j, k+1} = q_i \frac{A_a}{\Delta x \Delta y} = q_i \frac{(\Delta x - x)y}{\Delta x \Delta y} \end{eqnarray*}\]

\[ \rho_{j, k} = \sum_i q_i S({\bf R}_{j, k} - {\bf r}_i )\] \[ {\bf J}_{j, k} = \sum_i q_i {\bf v}_i S({\bf R}_{j, k} - {\bf r}_i )\]

2.18 直角坐标系插值

\[ {\bf F}_i(x, y) = \frac{1}{\Delta x \Delta y}(A_d {\bf F}_{j, k} + A_c {\bf F}_{j+1, k} + A_b {\bf F}_{j+1, k+1} + A_a {\bf F}_{j, k+1}) \]

2.19 直角坐标系插值

面积权重插值

面积权重插值

双线性插值

双线性插值

几何解释不同,结果一致。

2.20 \(r-\theta\)插值

\[ q_{j, k} = q_i \frac{(r^2_{j+1} - r_i^2)(\theta_{k+1} - \theta_i)}{(r^2_{j+1} - r_j^2)(\theta_{k+1} - \theta_k)} = q_i f_{j, k} \]

2.21 边界条件

边界条件

  • Periodic
  • Open
  • Reflecting

场:

  • 周期:\(F_{min}= F_{max}\)
  • 开放:PML
  • 反射:\(F = 0\)

3 epoch

3.1 epoch code flow

3.2 epoch结构

3.3 粒子运动 && 电流求解

3.4 epoch QED 物理过程

QED processes:

  • Radiation:
    • Synchrotron (nonlinear Compton scattering). \(e + n\omega \rightarrow e + \gamma\)

    • Bremsstrahlung \(e + Z \rightarrow e' + Z + \gamma\)

  • Pair production:
    • Breit-Wheeler pair production \(\gamma + n\omega \rightarrow e^- + e^+\)

    • Trident \(e + n\omega \rightarrow e + e^- + e^+\)

    • Bethe-Heitler pair production \(\gamma + Z \rightarrow e^- + e^+ + Z\)

3.5 QED algorithm

QED processes:

  • Synchrotron (nonlinear Compton scattering)
  • Breit-Wheeler pair production
  • process no needs for target

3.6 epoch-input

3.7 epoch-input

3.8 epoch-input

3.9 epoch-input

3.10 epoch-input

3.11 epoch-output

可以直接用Matlab脚本读出sdf文件中数据,输出为结构体形式:

data = GetDataSDF('filename.sdf')

之后,通过 data. 格式来访问数据中所有数据。

4 Hard coding

4.1 粒子结构

粒子数据结构

粒子数据结构一般为

typedef struct {
  x,y,z,px,py,pz,other
} ptclsdata;

粒子可以用平直数据表存储:

typedef struct {
int i;
ptclsdata myptcls[NMAX];
}

i代表当前存储的粒子数。

表的增加和删除
  • 增加粒子: \(i=i+1\),在表的末端添加一个粒子
  • 删除粒子:将被删除粒子和表的末端交换,\(i=i-1\)
  • 扩大表: 当\(i>NMAX*0.9\)
  • 减小表: 当\(i<NMAX*0.1\)

4.2 粒子数据

整体连续存储

在这种框架中,所有粒子被存储在一个大的表格中。优点是处理运动比较简单,缺点是在进行二体碰撞时需要 排序。即使不进行二体碰撞,也需要定期排序来提高计算速度。

按网格存储

每个网格维护一张表,共\(N_g\)个粒子表格,每当粒子运动越出网格边界时,就将其移动到正确的网格中。优点是 容易处理碰撞效应,而且无需排序即可达到较高的内存效率。缺点是粒子运动的计算量会增加。

并行处理

仅当使用区域分解算法时,在每个并行步结束时需要交换越界的粒子。

4.3 场量

Note
  • 全部场量 \(\rho\), \(j\), \(\phi\) 被全部进程所知道
  • 每个进程含有等量的粒子,推动粒子无须在进程中交换粒子信息
  • 每当推动粒子之后,电荷密度被广播到所有进程
  • 每次求解场方程,结果被广播到所有进程

模拟区域按照径向分割成子区域,每个子区域包含一定空间的场和粒子。当粒子移动越过边界时,将粒子数据传递到下一个节点。

4.4 场量

模拟区域按照径向分割成子区域,每个子区域包含一定空间的场和粒子。当粒子移动越过边界时,将粒子数据传递到下一个节点。

4.5 蒙特卡洛

PIC-MCC

在放电问题中,等离子体需要靠电离碰撞或者二次电子发射来进行维持。 为此,在PIC建模中必须引入描述带电粒子-中性粒子碰撞的过程。这可以通 过Monte-Carlo方法获得,这种用PIC方法描述粒子运动,用蒙特卡洛方法描述 粒子之间的近距离碰撞过程的方法称为PIC/MCC模拟。


抽样和蒙特卡洛

以碰撞电离过程为例,一个电子入射产生电离电子的几率为 \(n\sigma\delta t\),则对于一个宏粒子, 在 \(\delta t\) 时间内平均应该产生 \(nw\sigma\delta t\) 个次生电子,我们可以将其离散化为在 \(\delta t\) 时间内有 \(\nu=n\sigma\delta t\) 几率一个宏粒子变成两个宏粒子。

于是我们可以在每步推进中,增加一个蒙特卡洛抽样过程,即对每个宏粒子,产生随机数 \(r\) ,如果 \(r<\nu\) 则将 宏粒子变为两个,并随机设置两个粒子的出射能量动量(MC 碰撞)

4.6 蒙特卡洛

考虑一个相对于靶以速度\(V\)运动的入射粒子,粒子动能为E;靶粒子数密度为\(n\),而 入射粒子和靶之间碰撞的总截面为\(\sigma(E)=\sigma_1(E)+\sigma_2(E)+...\),下标代表各种不同类型的碰撞。那么,发生第i种碰撞的频率是 \(\nu_{coll}=n(x)\sigma_i V\),

在时间 \(\Delta t\)内,发生碰撞的总几率是 \[ \begin{equation} P=1-\exp[-V\Delta t \sigma(E)n(x)] \end{equation} \]

我们可以直接遍历所有粒子,并对每个粒子进行一次 Monte Carlo 抽样来完成模拟(无偏抽样),也可以 使用由 Vahedi 等人发明的空碰撞方法。

4.7 蒙特卡洛-空碰撞

为了提高性能,带电粒子和中性粒子之间的碰撞通常使用”空碰撞”方法描述:

首先设置最大碰撞几率,即 \[ \begin{eqnarray} P_{null}&=&1-\exp(-\nu_{max}\Delta t) \end{eqnarray} \]

设总的模拟粒子数为\(N\),计算出一个时间步内发生的最大可能碰撞数\(N_{coll}\)。 然后在总粒子中随机抽取\(N_{coll}\)个粒子,然后遍历抽取的样本粒子,对每个样本粒子,生成随机数\(R\),并根据\(r\)的具体大小执行如下操作: \[ \begin{eqnarray*} R&<&\frac {\nu_1(E)}{\nu_{max}} \qquad \text{Collision 1}\\ \frac {\nu_1(E)}{\nu_{max}} &<& R < \frac {\nu_1(E)+\nu_2(E)}{\nu_{max}}\ \ \ \ \ \ \text{Collision2} \\ &&....\\ \end{eqnarray*} \]

4.8 蒙特卡洛-末态

在碰撞中碰撞参数(瞄准距离)的设置是任意的,末速度利用出射角度计算 。由于散射几率存在一定的方向性,需要根据微分截面计算散射到各个方向的几率并按照这个几率抽样出射角度。在质心系中,出射到某一角度的几率是是 \[ \begin{equation} P(V,\chi)=\frac {\sigma(V,\chi)}{\sigma_{total}}\sin\chi d\chi d\phi \end{equation} \]

因此选择出射角度的方法是产生一个随机数\(r\),并且解方程 \[ \begin{eqnarray} \frac {2\pi}{\sigma_T}\int_0^{\chi}\sigma(V,\chi')\sin\chi'd\chi'=r \end{eqnarray} \] 设置出射角度后,出射粒子速度可以从质心系能量动量方程给出。

4.9 单体碰撞、多体碰撞

单体碰撞

带电粒子和背景气体的碰撞,包括弹性碰撞,激发和碰撞电离等,可以通过最简单的空碰撞抽样方法 完成,也可以实际遍历粒子进行抽样。在电离的情况下,设置出射粒子权重和入射粒子权重相等。

二体碰撞

二体碰撞包括电子-离子碰撞,电子-阳离子复合等等。主要的困难在于应该只允许同一网格内的粒子 发生二体碰撞,这可以通过设计Particle by Cell的数据结构完成,或者对粒子进行排序,或者遍历 所有粒子。

4.10 库伦碰撞

基本方法是所谓”binary collisions”:首先将将粒子按照网格分组,然后每个网格内的粒子随机配对,最后对于每一对粒子,进行一次随机碰撞。为了保证统计特性,碰撞出射角度分布要满足碰撞算子的结果: \[P(\chi)=\frac \chi {{\langle \chi^2 \rangle}_{\Delta t}}\exp\left(-\frac {\chi^2} {2{\langle \chi^2 \rangle}_{\Delta t}}\right)\] 其中出射角方差是 \[{\langle \chi^2 \rangle}_{\Delta t}=\frac {q_{\alpha}^2q_{\beta}^2}{2\pi\epsilon_0^2}\frac {n\Delta t \ln\Lambda_{\alpha\beta} }{\mu^2 V_{\alpha\beta}^3}\]

\(\mu\) 是折合质量,\(V_{\alpha\beta}\) 是相对速度,而 \(\ln \Lambda_{\alpha\beta}\) 是库仑对数,即 \(\Lambda_{\alpha\beta}= \frac {b_{max}}{b_{min}}\)\(b_{max}\)\(b_{min}\) 分别是最大和最小的瞄准距离(碰撞参数)。

4.11 库伦碰撞-Nanbu累积计算

Nanbu提出的累积碰撞模型稍稍降低了所需的计算量,思路是将多次小角度库仑碰撞归并为一次大角度的碰撞。由于这些次碰撞是随机的,可以按照随机过程累积为一次大角度碰撞。 于是不再需要限制 \(\Delta t\) 的长度。相反,\(\Delta t\) 用来计算库仑散射的分布,也就是散射几率 \(f(\chi)d\Omega\) 满足 \[ \begin{eqnarray} f(\chi)=\frac {A}{4\pi\sinh A}\exp(A \cos\chi)\\ \coth A-A^{-1}=\exp[-S_{\alpha\beta}(\Delta t)]\\ S_{\alpha\beta}(\Delta t)=\frac {\ln \Lambda_{\alpha\beta}}{4\pi}\left(\frac{q_{\alpha}q_{\beta}}{\epsilon_0\mu_{\alpha\beta}}\right)^2n_{\beta}V_{\alpha\beta}^{-3}\Delta t \end{eqnarray} \] 针对这个分布,Nanbu给出抽样公式为 \[ \begin{equation} \cos\chi=\frac 1A \ln [\exp(-A)+2R \sinh A] \end{equation} \]

4.12 库伦碰撞-朗之万近似

库仑碰撞的计算量实在太大,我们可以用郎之万近似来逼近它,也就是说,我们不再将粒子配对碰撞,而是直接将库仑碰撞作为噪声。 详细的内容见(W.M. Manheimer et al, J. Comput. Phys. 138 (1997) 563–584)

(近期进展:M. Sherlock, J. Comput. Phys. 227 (2008) 2286–2292,D.S. Lemons, D. Winske, W. Daughton, B. Albright, J. Comput. Phys. 228 (2009) 1391–1403, and A.M. Dimits et al. / Journal of Computational Physics 242 (2013) 561–580)

5 SLIPs-v1.0

5.1 SLIPs 框架

SLIPs框架

SLIPs框架

5.2 数据结构

SLIPs数据结构

SLIPs数据结构

5.3 框架

5.4 Run

while (t < domainInfo.tRange[1]) {
    if (simInfo.dumpInfo.use_dump) dumpData.dump(simBox, ...);
    if (simInfo.spinInfo.use_spin) spin->step(simBox, ...);
    for (auto &phys: sfqed) phys->step(simBox, ...);
    if (!biScatter.empty()) simBox.updateCellGroup();
    for (auto &phys: biScatter) phys->step(simBox, ...);
    pusher->step(simBox);
    if (!simInfo.pusherInfo.noCurrent && !domainInfo.noField)
        current->evaluate(simBox);//
    for (auto &l: laser) l.apply(simBox);
    solver->solve(simBox);
    simBox.istep++;
    t += domainInfo.dt;
}

5.5 粒子运动求解

\[\begin{eqnarray*} \mathbf{p}^{n-1/2} &=& \mathbf{p}^- - \frac{q\Delta t}{2m}\mathbf{E}^n, \label{boris1} \\ \mathbf{p}^{n+1/2} &=& \mathbf{p}^+ + \frac{q\Delta t}{2m}\mathbf{E}^n, \\ \mathbf{p}' &=& \mathbf{p}^- + \mathbf{p}^- \times \mathbf{t}, \\ \mathbf{p}^+ &=& \mathbf{p}^- + \mathbf{p}' \times \mathbf{s}, \\ \mathbf{t} &=& \frac{q\Delta t}{2m\gamma^n}\mathbf{B}^n, \\ \mathbf{s} &=& \frac{2\mathbf{t}}{1+t^2}, \label{boris6} \end{eqnarray*}\]

where \(\gamma^n = \sqrt{1+\mathbf{p}^2}\).

5.6 粒子运动求解器

在SLIPs中粒子运动求解器初始化 simulator.cpp => init():


if (pusherInfo.useImplicit) 
    pusher = new(implicitLorentz);
else {
    if (pusherInfo.useBoris)
        pusher = new (_Boris);
    else {
        if (pusherInfo.useLandauLifshitzPush) {
            pusher = new (_Landau);
        } else if (pusherInfo.useModifiedLLPush) {
            pusher = new (_MLandau);
        } else pusher = nullptr;
    }
}
if (pusher != nullptr) pusher->init(simBox, ...);

5.7 粒子运动求解器

在pusher是一个粒子推进基类(虚),真正的粒子推进必须由该类继承,并实现:

#ifndef SLIPS_BORIS_H
#define SLIPS_BORIS_H
#include "pusher.h"
namespace slips {
    class _Boris : public __Pusher {
    public:
        void push(_particle &particle, _simBox &simBox);
        arma::vec3 uminus, t, uprime, s, uplus, ufinal;
    };
} // slips
#endif //SLIPS_BORIS_H

5.8 粒子运动求解器

  • 在pusher是一个粒子推进基类(虚),真正的粒子推进必须由该类继承,并实现。 例如,目前已实现的包含 - 无辐射修正的Newton-Lorentz为Boris: \[\frac{d{\mathbf{p}}}{d t} = \frac{q}{m}(\mathbf{E} + \beta \times \mathbf{B})\]

  • 经典辐射修正的Newton-Landau-Lifshitz为LL: \[\begin{equation} \begin{split} \mathbf{F}_{\mathrm{RR, classical}}& =\frac{2e^3}{3mc^3} \bigg \lbrace \gamma \bigg[ \bigg (\frac{\partial}{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot\nabla \bigg )\mathbf{E}+\frac{\mathbf{p}}{\gamma mc} \times \bigg (\frac{\partial}{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot\nabla \bigg )\mathbf{B} \bigg] \\ & + \frac{e}{mc} \bigg[\mathbf{E}\times \mathbf{B}+\frac{1}{\gamma mc}\mathbf{B}\times (\mathbf{B}\times \mathbf{p})+\frac{1}{\gamma mc}\mathbf{E}(\mathbf{p \cdot E}) \bigg] \\ & -\frac{e\gamma}{m^2c^2}\mathbf{p} \bigg [ \bigg(\mathbf{E}+\frac{\mathbf{p}}{\gamma mc}\times \mathbf{B} \bigg)^2-\frac{1}{\gamma^2m^2c^2}(\mathbf{E \cdot p})^2 \bigg ] \bigg\rbrace \end{split} \end{equation} \]

5.9 粒子运动求解器

量子辐射修正的Newton-MLL为: \(\mathbf{F}_\mathrm{QRR}=q(\chi)\mathbf{F}_\mathrm{RR}\),其中

\[\begin{equation*} q(\chi_e) \approx \left[1+4.8(1+\chi_e)\ln(1+1.7\chi_e)+2.44\chi_e^2\right]^{-2/3}, \end{equation*}\]
\[\begin{equation*} \mathbf{p}^{n+1/2} = \mathbf{p}^{n+1/2}_L + \mathbf{p}^{n+1/2}_R - \mathbf{p}^{n-1/2} = \mathbf{p}^{n+1/2}_L + \mathbf{f}_R^n \Delta t. \end{equation*}\]
\[\begin{eqnarray*} q(\chi) &=& \frac{I_\mathrm{QED}}{I_\mathrm{C}}, \\ I_\mathrm{C} &=& \frac{2e^4E'^2}{3m^2c^3}, \\ I_\mathrm{QED} &=& mc^2\int c(k\cdot k')\frac{dW_{fi}}{d\eta d r_0}dr_0. \end{eqnarray*}\]

5.10 粒子运动求解器

Boris pusher的实现:

Boris pusher

Boris pusher

5.11 粒子运动求解器

Landau-Lifshitz pusher的实现:

LL pusher

LL pusher

5.12 求解器

pusher

pusher

solver

solver

5.13 电磁场求解

3D case:

\[\begin{eqnarray} \begin{split} \frac{E_y^{n+1}-E_y^n}{\Delta t}\bigg|_{i+1/2, j, k+1/2} = & -\frac{B_{i+1, j, k + 1/2} - B_{i, j, k + 1/2}}{\Delta x}\bigg|_z^{n+1/2} + \\ & \frac{B_{i+1/2,j, k+1} - B_{i+1/2,j, k}}{\Delta z}\bigg|_x^{n+1/2}- \\ & J_{y,i+1/2,j, k+1/2}^{n+1/2} \nonumber \\ \frac{B_z^{n+1/2}-B_z^{n-1/2}}{\Delta t}\bigg|_{i, j, k+1/2} = & \frac{E_{i+1/2, j, k} - E_{i-1/2, j, k}}{\Delta x}\bigg|_y^{n} - \\ & \frac{E_{i+1/2,j+1, k} - E_{i+1/2,j, k}}{\Delta y}\bigg|_x^{n} \nonumber \end{split}, \end{eqnarray}\]

5.14 如何引入新的电磁场求解

Yee 二阶

5.15 如何引入新的电磁场求解

初始化

5.16 非线性康普顿散射

5.17 非线性Breit-Wheeler散射

5.18 如何引入两体散射过程

5.19 如何引入两体散射过程

5.20 如何引入两体散射过程

5.21 输入文件

SLIPs/i-SLIPs的输入文件采用toml格式,

5.22 输入文件

[domainInfo]
ncell_global = [512, 1, 1]
xyzmin_global = [-25.13, 0.0, 0.0]
xyzmax_global = [25.13, 1.0, 1.0]
ngrid_ghost = [4, 0, 0]
tRange = [0.0, 100.0]
dt = 0.02
printStep = 100
useProcessBar = false
noField = false

[[laserInfo]]
useLaser = false
boundary = "XMin"
ex = "0"
bx = "0"
ez = "0"
by = "0"
ey = "
  exp(-((t - 23.0) / 10.0)^4 
  - (y/ 5.0)^2) * 5.0"
bz = "
  exp(-((t - 23.0) / 10.0)^4 
  - (y/ 5.0)^2) * 5.0"

5.23 输入文件

5.24 场信息

5.25 配置文件读入

5.26 读入配置文件

5.27 配置文件读入

5.28 编译

5.29 编译 & 运行

$ cmake path_to_src_top
$ make -jN
$ make install
mpirun -np N path_to_bin/SLIPs input.toml

数据合并(utils):

python ../utils/mergeField.py path_to_h5 prefix

5.30 下一步?

5.31 Thanks

Back to top