Implicit discretization schemes for full-kinetic ion and drift-kinetic electron simulationsJournal: Journal of Computational Physics

In this subsection, we define the notation used throughout the paper for both field and particle quantities. A summary is provided in Table 1 for quick reference.

The FIDES model employs a full-kinetic description for ions, governed by the Vlasov equation

while the electrons are described by the drift kinetic equation

The guiding center velocity is

with $\mathbf{b}=\mathbf{B}_{0}/B_{0}$. Here $\mathbf{v}_{D}$ contains the grad-B and curvature drift, and $\mathbf{v}_{E}$ is the $\mathbf{E}\times\mathbf{B}$ drift,

where the subscript 1 denotes perturbed quantities. The perturbed magnetic field is assumed to be much smaller than the equilibrium magnetic field. Therefore, the amplitude of the total magnetic field is $|\mathbf{B}_{0}+\mathbf{B}_{1}|\approx B_{0}+B_{1\|}$ and the direction of the total magnetic field is $\left(\mathbf{B}_{0}+\mathbf{B}_{1}\right)/|\mathbf{B}_{0}+\mathbf{B}_{1}|\approx\mathbf{b}+\mathbf{B}_{1\perp}/B_{0}$ in Eq. (3). The time derivative of the electron kinetic energy $\varepsilon_{e}=m_{e}v^{2}/2$ is

where the magnetic moment $\mu_{e}=m_{e}v_{\perp}^{2}/2B_{0}$ is an adiabatic invariant in the drift kinetic model.

To reduce the particle noise in simulations, we employ the $\delta f$ method for both species in the present work [12]. However, $\delta f$ scheme suffers from severe numerical problems when the perturbed distribution function is comparable to the equilibrium distribution function. Simulating large perturbation problems lies beyond the scope of the present work.

Regarding the field equations, the implicit discretization scheme is built upon a reduced form of the Maxwell equations

where the displacement current term in Ampere’s law is dropped. This omission is motivated by the need to numerically eliminate plasma oscillations at the plasma frequency $\omega_{pe}=\sqrt{n_{e}e^{2}/m_{e}\varepsilon_{0}}$ (with $\varepsilon_{0}$ the vacuum permittivity), which can be comparable to the electron cyclotron frequency. The plasma oscillations are not relevant to the physics of interest and, if retained, would violate the ordering assumptions of the electron drift kinetic equation.

An important consequence of this reduced form of Ampere’s law is that it inherently ensures the quasi-neutrality condition. Taking the divergence of Ampere’s law and noting that the divergence of a curl is zero, we obtain $\nabla\cdot\mathbf{J}=0$. Meanwhile, the charge conservation equation holds as $\partial(q_{i}n_{i}-en_{e})/\partial t+\nabla\cdot{\mathbf{J}}=0$, which can be derived from the first-order moment of the ion Vlasov equation and electron drift kinetic equation. Combining $\nabla\cdot\mathbf{J}=0$ with the charge conservation equation yields $\partial(q_{i}n_{i}-en_{e})/\partial t=0$. Thus, within this formulation, quasi-neutrality is consistently maintained.

The original Maxwell equations (6) are ill-posed for explicit particle-pushing schemes. This is attributed to the fact that particle noise breaks the $\nabla\cdot\mathbf{J}=0$ condition required by Ampere’s law, resulting in inaccurate electromagnetic fields. A common solution is to transform the Ampere’s law into the generalized Ohm’s law [8, 9, 10, 11]. Here, in the context of particle-in-cell (PIC) magnetic confinement fusion (MCF) simulations, the generalized Ohm’s law is used to solve for the electric field using the current computed from particles. This usage differs from the classical Ohm’s law, which serves as a constitutive relation expressing the current in terms of the electric field. Depending on the direction relative to the magnetic field, the generalized Ohm’s law can be split into a parallel component and a perpendicular component. In previous work [8, 9, 10, 11], both components have been used as field equations. However, a numerical cancellation problem arises in the parallel Ohm’s law when electrons behave adiabatically for $k_{\|}v_{te}\gg\omega$ [8, 10, 11]. In this regime, the leading-order terms $E_{1\|}$ and $-\nabla_{\|}p_{1,e\|}$ nearly cancel, and governing dynamics appears at the next order. Consequently, numerical noise in $p_{1,e\|}$ can lead to a severe loss of accuracy.

To address the numerical cancellation problem inherent in the parallel Ohm’s law, we propose using the Ampere’s law directly to solve for the electric field. Since $E_{\|}$ does not appear explicitly in Ampere’s law but enters implicitly through the electron response, an implicit treatment of the $E_{\|}$ contribution to the electron response has to be used. This leads to the formulation of an implicit parallel Ampere’s law, which employs a lower-order electron velocity moment and, as shown in Section 3, effectively mitigates the cancellation problem. In our algorithm, we adopt this approach by replacing the parallel Ohm’s law with the implicit parallel Ampere’s law, while retaining the perpendicular Ohm’s law. The electric field is then obtained by solving these two equations simultaneously as a coupled system.

To develop an implicit scheme for the parallel Ampere’s law, we advance electrons using an implicit $E_{\|}$ scheme

which have been normalized according to the conventions in section 2.1.3. Here, the intermediate electron weight $w_{e}^{*}$ is first advanced from $w_{e}^{n}$ via the first equation. Once the field equations are solved, the electron weight $w_{e}^{n+1}$ at time step $t^{n+1}$ is obtained by advancing $w_{e}^{*}$ using the updated parallel electric field $E_{1\|}^{n+1}$. Since the magnetic moment $\mu_{e}$ is a constant of motion, its superscript $n$ has been omitted in Eq. (7). Under the implicit $E_{\|}$ scheme, the electron parallel current $J_{e\|}^{n+1}$ contains two parts

where $J_{e\|}^{*}$ and the second term of $J_{e\|}^{n+1}$ are, respectively, the intermediate and implicit parts in the electron parallel current. In Eq. (8), $S\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)$ is the shape function, where $\mathbf{x}_{g}$ and $\mathbf{x}_{ej}^{n+1}$ are the positions of the grid point and the electron, respectively. In the FIDES algorithm, the particle positions and velocities at $t^{n+1}$ are computed before accumulating particle moments. In linear simulations, the particle trajectories are advanced along unperturbed orbits, independently of the perturbed fields. In nonlinear simulations, the particle positions and velocities are updated iteratively using the fields from previous iteration, so that $v_{ej\|}^{n+1}$ and $\mathbf{x}_{ej}^{n+1}$ are available when the particle moments are computed.

Substituting the electron parallel current (8) into the parallel Ampere’s law of Eq. (6), we can get

where $\beta_{e}=\mu_{0}n_{e}T_{e}/B_{0}^{2}$ and $N_{p}$ is the particle number in one grid. Next, substituting $\mathbf{B}_{1}^{n+1}$ with $\mathbf{B}_{1}^{n}-\Delta t\nabla\times\mathbf{E}_{1}^{n+1}$ and collecting all terms that depend on the electric field at time step $t^{n+1}$ on the left-hand side, we obtain the implicit parallel Ampere’s law, which serves as one of the equations for the electric field,

After computing the ion and electron currents, and given $\mathbf{B}_{1}^{n}$, this equation is solved together with the perpendicular Ohm’s law (introduced later) to obtain the electric field at time step $t^{n+1}$. Note that regardless of whether the field equations are solved using a finite difference method or a spectral method, the left-hand side of Eq. (10) contains a discrete summation term over electrons, $\left(\beta_{e}/N_{p}\right)\sum_{j}\Delta t\left(v_{ej\|}^{n+1}\right)^{2}E_{1\|}^{n+1}\left(\mathbf{x}_{ej}^{n+1}\right)S\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)$. This term renders the discretized coefficient matrix of field equations complex and time-step-dependent, precluding the use of cost-saving techniques such as performing a single LU decomposition and storing the matrix for use in later time steps. Therefore, in FIDES, we propose to solve the implicit parallel Ampere’s law and perpendicular Ohm’s law in an iterative manner. In particular, we first find an approximate expression for $\left(\beta_{e}/N_{p}\right)\sum_{j}\Delta t\left(v_{ej\|}^{n+1}\right)^{2}E_{1\|}^{n+1}\left(\mathbf{x}_{ej}^{n+1}\right)S\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)$,

Here, we have approximated the shape function as $S\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)\approx\Delta V\delta\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)$, where $\Delta V$ is the volume of a grid cell. This approximation becomes increasingly accurate as the grid size decreases. In the limit of a sufficiently large number of marker particles, we can approximate $\left(\Delta V/N_{p}\right)\sum_{j}\delta\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)\delta\left(v_{\|}-v_{ej\|}^{n+1}\right)$ as the normalized marker distribution $exp[-m_{e}v_{\|}^{2}/\left(2m_{i}\right)]/(2\pi m_{i}/m_{e})^{0.5}$ in Eq. (11),

Substituting Eq. (12) into Eq. (10), treating the difference between $\left(\beta_{e}/N_{p}\right)\sum_{j}\Delta t\left(v_{ej\|}^{n+1}\right)^{2}E_{1\|}^{n+1}\left(\mathbf{x}_{ej}^{n+1}\right)S\left(\mathbf{x}_{g}-\mathbf{x}_{ej}^{n+1}\right)$ and $\Delta t\beta_{e}\left({m_{i}}/{m_{e}}\right)E_{1\|}^{n+1}\left(\mathbf{x}_{g}\right)$ as a small perturbation and moving them to the right-hand side, we can get the iterative form of the implicit parallel Ampere’s law,

where we have updated the superscript notation to distinguish between time step and iteration: the electric field at time step $n+1$ is now denoted with iteration indices $k$ and $k+1$ within the iterative loop. Specifically, for the initial guess $(k=0)$, we set $\mathbf{E}_{1}^{k}=\mathbf{E}_{1}^{n}$. In FIDES, the iteration loop terminates when the relative change between successive iterates falls below a prescribed threshold, $\|\mathbf{E}_{1d}^{k+1}-\mathbf{E}_{1d}^{k}\|_{2}/\|\mathbf{E}_{1d}^{k}\|_{2}<\text{tol}$, $(d=x,~y,~z)$, where $||\cdot||_{2}$ denotes the global 2-norm (i.e., the Euclidean norm of the vector formed by collecting the electric field values at all grid points). The converged solution is then assigned as $\mathbf{E}_{1}^{n+1}=\mathbf{E}_{1}^{k+1}$. When the particle number and grid number are sufficient, convergence is typically achieved within a few steps; a convergence test can be seen in Fig. 3. It is important to note that for quantities which remain unchanged during the iteration loop, such as $\mathbf{B}_{1}^{n}$ in Eq. (13), the superscript continues to denote the time step n, not the iteration index.

Upon convergence, the iterative form of the implicit parallel Ampere’s law, Eq. (13), reverts to its original formulation, Eq. (10). Comparing the two equations reveals that the coefficient matrix in Eq. (13) is considerably simpler and, importantly, independent of the time step. This allows us to perform a single LU decomposition and store the matrix, significantly reducing the computational cost of solving the field equations.

For typical parameters in fusion plasmas, the presence of high-frequency waves with perpendicular electric fields (e.g., compressional Alfvén, ion Bernstein, and extraordinary waves) forces an extremely small timestep $\Omega_{ci}\Delta t<0.01$. To eliminate this strict limitation, we push ions in an implicit $\mathbf{E}_{\perp}$ scheme [8, 9]

Here, the ion equilibrium distribution is assumed as a local Maxwellian distribution and we can get the ion perpendicular current

Given that the perturbed electron distribution is independent of gyrophase in the drift kinetic model, it follows that the perpendicular electron current is alternatively determined from the electron perpendicular momentum equation

where $p_{1,e\perp}^{n+1}$ is the perturbed electron perpendicular pressure.

Substituting Eq. (15) and (16) into the perpendicular Ampere’s law of Eq. (6), we can get

Then, taking the cross product of this equation with $\mathbf{b}$, substituting the magnetic field $\mathbf{B}_{1}^{n+1}$ using Faraday’s law, and moving all terms involving $\mathbf{E}_{1}^{n+1}$ to the left-hand side, we obtain the implicit perpendicular Ohm’s law.

For reasons similar to those discussed for the implicit parallel Ampere’s law, Eq. (18) is also solved iteratively in FIDES. Following the same approach of Eq. (11), we introduce the approximation $S\left(\mathbf{x}_{g}-\mathbf{x}_{ij}^{n+1}\right)\approx\Delta V\delta\left(\mathbf{x}_{g}-\mathbf{x}_{ij}^{n+1}\right)$. Under the assumption of a sufficiently large number of marker particles, $\left(\Delta V/N_{p}\right)\sum_{j}\delta\left(\mathbf{x}_{g}-\mathbf{x}_{ij}^{n+1}\right)\delta\left(\mathbf{v}-\mathbf{v}_{ij}^{n+1}\right)$ is then approximated by the normalized ion marker distribution $exp\left[-(T_{e}v^{2})/(2T_{i})\right]/(2\pi T_{i}/T_{e})^{1.5}$, yielding

It should be mentioned that the perturbed electron perpendicular pressure $p_{1,e\perp}^{n+1}$ also contains an implicit contribution, arising from the implicit $E_{\|}$ scheme used to advance electron weights (see Eq. (7)),

Following a derivation analogous to that of Eq. (19), the second term (the implicit part) on the right-hand side of Eq. (20) can be approximated by zero. This approximation does not affect the iterative form of Eq. (18). Therefore, we have not expanded the expression for $p_{1,e\perp}^{n+1}$ in Eq. (18). However, in the actual FIDES implementation, the implicit part of $p_{1,e\perp}^{n+1}$ is retained and computed iteratively. This is because the particle weights are updated after each field equation iteration, and the particle moments including $p_{1,e\perp}^{n+1}$ are subsequently recalculated as a whole.

Substituting Eq. (19) into Eq. (18) and moving the difference between the discrete summation term and its approximation term to the right-hand side, we obtain the iterative form of the implicit perpendicular Ohm’s law,

where the superscript of the the electric field has been changed from time step $n+1$ to iteration indices $k$ and $k+1$. The iteration process of Eq. (21) is similar to that of Eq. (13).

In FIDES, the implicit parallel Ampere’s law (13) and the implicit perpendicular Ohm’s law (21) are solved simultaneously as a coupled system to obtain the electric field. The numerical method employed depends on the complexity of the equilibrium profile. For problems with a uniform equilibrium, a spectral method is used. For cases where equilibrium non-uniformities are restricted to the x direction, we apply a hybrid approach that combines a finite difference discretization in x with a Fourier decomposition in the y and z directions. For more complex problems, the dimension of the coefficient matrix becomes large. In such cases, cost-efficient methods are required, such as precomputing and storing an LU decomposition or employing a Krylov subspace method. In this paper, for the purpose of demonstration, we present the spectral method for solving Eq. (13) and Eq. (21). After Fourier transformation in x, y and z directions, the field equations become

where elements of the coefficient matrix are

On the right-hand side of Eq. (22), $\tilde{r}_{x}\left(\mathbf{k}\right),\tilde{r}_{y}\left(\mathbf{k}\right)$ and $\tilde{r}_{z}\left(\mathbf{k}\right)$ can be expressed as

where $\tilde{\left(\cdot\right)}$ and $F\left\{\cdot\right\}$ denote Fourier transformed quantities and $\mathbf{k}=k_{x}\hat{\mathbf{x}}+k_{y}\hat{\mathbf{y}}+k_{z}\hat{\mathbf{z}}$ is the wavenumber. At each field iteration, the updated electric field $\tilde{\mathbf{E}}_{1}^{k+1}\left(\mathbf{k}\right)$ is obtained by solving Eq. (22). This Fourier-space solution is then filtered and transformed back to real space via an inverse Fourier transform, yielding $\mathbf{E}_{1}^{k+1}\left(\mathbf{x}_{g}\right)$ at the grid points.