A conservative implicit-PIC scheme for the hybrid kinetic-ion fluid-electron plasma model on curvilinear meshes

For the time advance, we employ the second-order implicit-midpoint method. For a variable $\chi$ with known value at timestep $t=n\Delta t$, i.e. $\chi^{n}$, the new value $\chi^{n+1}$ is found implicitly from $(\chi^{n+1}-\chi^{n})/\Delta t=f(\chi^{n+1/2})$ where $\chi^{n+1/2}=\tfrac{1}{2}(\chi^{n+1}+\chi^{n})$. We note that this differs from the Crank-Nicolson method in the discretization of non-linear terms. This timestepping scheme is used here both for the grid-based equations (Sec. 3.3), and for the particle advance (Sec. 3.4) where it gives long term accuracy for ion orbit integration (it is symplectic). This method was found to give exact discrete conservation in the Cartesian formulation of this hybrid-PIC kinetic-ion and fluid-electron scheme [1].

Assuming a static and non-singular curvilinear map from logical space $\xi^{\alpha}=(\xi,\eta,\mu)$ to physical space $\boldsymbol{x}=(x,y,z)$, i.e. $\boldsymbol{x}(\xi^{\alpha})$, the Jacobian of the mapping is defined as $J=|\partial\boldsymbol{x}/\partial\xi^{\alpha}|$, the contravariant basis vectors as $\boldsymbol{\nabla}\xi^{\alpha}$, and the covariant basis vectors as $\partial_{\alpha}\boldsymbol{x}\equiv\partial\boldsymbol{x}/\partial\xi^{\alpha}$. Then, for a generic vector $\boldsymbol{S}$, the contravariant components are defined as $S^{\alpha}=\boldsymbol{S}\cdot\boldsymbol{\nabla}\xi^{\alpha}$, and the covariant components as $S_{\alpha}=\boldsymbol{S}\cdot\partial_{\alpha}\boldsymbol{x}$. Indices can be lowered as $S_{\alpha}=g_{\alpha\beta}S^{\beta}$ via the metric tensor $g_{\alpha\beta}\equiv\partial_{\alpha}\boldsymbol{x}\cdot\partial_{\beta}\boldsymbol{x}$, and raised as $S^{\alpha}=g^{\alpha\beta}S_{\beta}$ where $g_{\alpha\beta}g^{\beta\gamma}=\delta_{\alpha}^{\gamma}$. A gives some further definitions and vector identities used herein.

In hybrid-PIC discretizations of the model described in Section 2, Eqs. (2, 3, 4) for $\boldsymbol{B}$, $\boldsymbol{E}$, and $p_{e}$ are typically solved on a spatial grid. Here, we use a uniform structured grid in logical space with all grid variables $\chi_{g}\equiv\chi_{ijk}$ defined at the position of the cell centers in logical space $(\xi_{i},\eta_{j},\mu_{k})$ for $i\in[1,N_{\xi}]$, $j\in[1,N_{\eta}]$, $k\in[1,N_{\mu}]$. The subscript $g$ (shorthand for the 3-dimensional grid index $ijk$) is used to denote grid-based variables herein. The spatial derivatives of these mesh variables are computed using a second-order centered finite-difference approximation, as e.g. $(\partial_{\xi}\chi)_{g}\equiv(\partial_{\xi}\chi)_{ijk}=(\chi_{i+1jk}-\chi_{i-1jk})/(2\Delta\xi)$, where $\Delta\xi$ is the uniform grid spacing on the logical mesh.

The covariant electric field $E_{\alpha,g}^{n+1/2}$ is found via static evaluation of Ohm’s law at the $n+1/2$ time level as

where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol. The density and velocity moments of the kinetic ions, $n_{g}$ and $(u^{\beta})_{g}$, are defined below. $(B^{\gamma})_{g}$ are the contravariant components of the magnetic field, $(j^{\beta})_{g}=(1/J_{g})\epsilon^{\beta\sigma\alpha}(\partial_{\sigma}(g_{\alpha\gamma}B^{\gamma}))_{g}$ are the contravariant components of the current density, and $(p_{e})_{g}$ is the electron scalar pressure.

The contravariant magnetic field components are calculated via Faraday’s equation as

The scalar electron pressure is updated as

where $(u_{e}^{\alpha})_{g}=(u^{\alpha})_{g}-(j^{\alpha})_{g}/en_{g}$ is the electron bulk velocity. Note that we have multiplied Eq. (8) by $J_{g}$, which we assume to be constant in time. We note that these choices for discretization of the spatial grid equations are motivated by those used for a non-staggered resistive-MHD scheme described in Ref. [44].

To close the set of equations, electromagnetic fields are scattered from the spatial grid to the particle positions, and moments of the particle distributions are gathered for use in Eqs. (6-8). For the momentum conservation properties of our scheme described in Sections 4.2.1-4.2.2, we find it necessary to scatter the Cartesian components of the electromagnetic fields to the particle positions as

where $S(\xi^{\alpha}_{g}-\xi^{\alpha}_{p})\equiv S_{m}(\xi_{i}-\xi_{p})S_{m}(\eta_{j}-\eta_{p})S_{m}(\mu_{k}-\mu_{p})$ is the tensor product of $m$-th order B-spline shape functions. The electromagnetic fields are assumed to be constant in time over the particle substeps (so are given superscript $n+1/2$), such that the change in force occurs only due to the change in particle substep position. The Cartesian vector values of the fields are found on the spatial mesh as $\boldsymbol{E}_{g}=E_{\alpha,g}\left(\boldsymbol{\nabla}\xi^{\alpha}\right)_{g}$, $\boldsymbol{B}_{g}=(B^{\alpha})_{g}\left(\frac{\partial\boldsymbol{x}}{\partial\xi^{\alpha}}\right)_{g}$ prior to interpolation.

The moments are gathered using the orbit-averaging technique [46], where a contribution to the moments at each sub-step is weighted by the fraction of the macro-scale timestep, $\Delta\tau_{p}^{\nu}/\Delta t$. This orbit averaging technique increases the number of statistical samples in the calculation of the moments, thus reducing the noise. Also, as shown below, it is required to give discrete conservation properties in the presence of subcycling. The quasi-neutral density is gathered as

where $w_{p}$ is the macro-particle weight, which we assume to be constant in time for each particle, and $\Delta V=\Delta\xi\Delta\eta\Delta\mu$ is the constant logical cell volume. The Cartesian bulk momentum is gathered as

The contravariant component used in the field and moment equations is then found as $(nu^{\alpha})_{g}^{n+1/2}=(n\boldsymbol{u})_{g}\cdot(\boldsymbol{\nabla}\xi^{\alpha})_{g}$.

The discrete linear momentum is defined as $\boldsymbol{P}^{n+1}=\sum_{p}w_{p}m_{p}\boldsymbol{v}_{p}^{n+1}$. The rate of change over a macro-scale timestep can then be written as

Using Eqs. (11,  12, 13, 14 and 15), this can be written as

which can be written in curvilinear form as

Then, using the covariant form of Ohms law from Eq. (6), and cancelling the convective electric field terms gives

Finally, using the definition of the current density $\left(j^{\beta}\right)_{g}^{n+1/2}$ and the contraction of indices for the product of Levi-Cevita symbols $\epsilon_{\alpha\beta\gamma}\epsilon^{\beta\sigma\kappa}=-(\delta_{\alpha}^{\sigma}\delta_{\gamma}^{\kappa}-\delta_{\alpha}^{\kappa}\delta_{\gamma}^{\sigma})$ gives

As in Ref. [1], our implementation solves for the magnetic vector potential $(A_{\alpha})^{n+1}_{g}$ rather than the magnetic field formulation $(B^{\alpha})_{g}^{n+1}$ used in Eq. (7). The Weyl gauge is chosen with electrostatic potential $\phi=0$, such that $(A_{\alpha})^{n+1}_{g}$ is found from the solution of

We note that all of the conservation properties of Sec. 4 also hold for the $\boldsymbol{A}$-formulation with the definition of $(B^{\alpha})_{g}=\epsilon^{\alpha\beta\gamma}\left(\partial_{\beta}A_{\gamma}\right)_{g}/J_{g}$.

The non-linear algebraic equations (43, 8) can be written in the form $\boldsymbol{G}(\boldsymbol{y}^{n+1})=\boldsymbol{0}$ for the solution vector $\boldsymbol{y}^{n+1}=\left((A_{\alpha})_{g}^{n+1},(p_{e})_{g}^{n+1}\right)$. The full details of the method of solution closely follow Ref. [1]. The system is approximately inverted by using a Jacobian-free Newton Krylov (JFNK) method [49] using flexible-GMRES [50] as a linear solver. The outer Newton iteration is iterated until the residual is converged to a user specified non-linear tolerance $\epsilon_{t}$. Note that the particle quantities $(\xi^{\alpha}_{p},\boldsymbol{v}_{p})$ do not enter the residual directly; rather, the particles are advanced by Eqs. (9, 11) and moments gathered by Eqs. (14, 15) in each residual evaluation. To do this, a Picard-iterated implicit Boris method is used to solve Eqs. (9, 11) - see Ref. [1].

Effective preconditioning plays a primary role in the efficiency of methods such as these. However, the discussion of a suitable physics-based preconditioner is outside of the scope of the present paper, and will be documented elsewhere. For these results, we do not perform preconditioning and, as such, we typically use a small timestep and do not consider performance data for the results presented here.

The discrete model is implemented in normalized form using natural magnetized ion units. For a typical magnetic field strength $B_{0}$, density $n_{0}$, and ion mass $m_{0}$, velocities are normalized by the Alfvén speed $v_{0}=v_{A0}=B_{0}/\sqrt{m_{0}\mu_{0}n_{0}}$, times by the inverse cyclotron period $t_{0}=\Omega_{c0}^{-1}=m_{0}/eB_{0}$, and lengths by the ion skin depth $L_{0}=d_{i0}=v_{A0}/\Omega_{c0}$. The different species of kinetic ions are distinguished by their normalized charge $Z_{s}=q_{s}/e$, mass $M_{s}=m_{s}/m_{0}$, and density $N_{s}=n_{s0}/n_{0}$. If they are initialized with a Maxwellian distribution function with temperature $T_{s0}$, we define the thermal velocity as $v_{Ts0}=\sqrt{T_{s0}/m_{s}}$ and plasma-beta $\beta_{s0}=2M_{s}N_{s}(v_{Ts}/v_{A0})^{2}=2\mu_{0}n_{s0}T_{s0}/B_{0}^{2}$. The electron sound speed is defined as $C_{s}=\sqrt{T_{e0}/m_{0}}$ and the ratio of ion to electron temperatures is $\tau_{s}=T_{s0}/T_{e0}$.

In the first numerical example, we consider the ion Landau damping of an ion acoustic wave using a 2D non-orthogonal curvilinear mesh. This is an electrostatic problem and the motivation is to numerically verify the novel conservation properties of the scheme for this limit (momentum and energy conservation, see Sec. 4.2.1) for non-orthogonal mapped meshes. The physical problem set-up considered is similar to that described in Ref. [1]. The simulation is initialized with a single ion species with $Z_{i}=M_{i}=1$ with a Maxwellian velocity function with thermal speed $v_{Ti0}=\sqrt{1/3}$. We set $\gamma=5/3$ and use a temperature ratio of $\tau=0.2$. The ions are given a perturbation to their velocity such that the initial bulk velocity is given by $\boldsymbol{\delta u}_{i}=0.01\cos{\left(k_{x}x+k_{y}y\right)}\left(\boldsymbol{\hat{x}}+\boldsymbol{\hat{y}}\right)$ with $k_{x}=k_{y}=\pi/8$. We use a non-orthogonal sinusoidal mesh with periodic boundary conditions for this problem, specified by the map

Here, the parameter $\sigma$ describes the distortion of the mesh. Fig. 1 shows the initial conditions for the $\boldsymbol{\hat{x}}$-component of the ion velocity moment, along with the spatial mesh used for $\sigma=1$.

Additional numerical parameters used are $64\times 64$ cells, with $5\times 10^{4}$ particles/cell initialized with a quasi-quiet start [51] from a low discrepancy Hammersley sequence [52] to reduce noise. We use second-order quadratic-spline shape functions and apply two passes of conservative binomial smoothing, see C. The timestep used is $\Delta t=0.02$, and the non-linear tolerance used is $\epsilon_{t}=10^{-12}$.

Figure 2 (left) shows the amplitude of the ion velocity moment vs time for the sinusoidal mesh (red) of Eqs. (44, 45), as well as a uniform Cartesian mesh with the same resolution (blue). The amplitude of the wave decreases exponentially with time due to kinetic ion Landau damping. The theoretical damping rate is found from the dispersion relation $Z^{\prime}(\xi)=2\tau$, where $Z$ is the plasma dispersion function and $\xi=(\omega+i\gamma)/k\sqrt{2}v_{Ti0}$ is the normalized complex frequency. For our initial conditions, the theoretical value is $\gamma=-0.0744636$, which is shown as the black dashed line on Fig. 2 (left). There is good agreement for both the Cartesian and curvilinear mesh cases with this theoretical damping rate until the noise due to a finite number of particles overcomes the signal.

The conservation errors are shown in Fig. 2 (right) for both the Cartesian mesh (blue) and curvilinear mesh (red) simulations. The momenta and energy are conserved to better than $10^{-13}$ and $10^{-12}$, respectively, which is on the level of the non-linear tolerance ($10^{-12}$) used. This provides numerical verification of the result of Section 4.2.1, that, for the electrostatic limit, the total linear momentum and energy is conserved for general curvilinear mapped meshes.

In the second numerical example, we simulate the right-hand polarized Alfvén-whistler wave in a cold and uniform background plasma. This wave has an exact analytical solution to the cold-plasma hybrid model for arbitrary perturbation amplitude, and therefore can be used to perform a formal grid convergence study. All simulations are performed in 2D, with the initial background plasma specified by $Z_{i}=M_{i}=N_{i}=1$, $T_{i}=T_{e}=0$, and the magnetic field oriented diagonally in physical space as $\boldsymbol{B}/B_{0}=(\boldsymbol{\hat{x}}+\boldsymbol{\hat{y}})/\sqrt{2}$.

The eigenmode of the wave is perturbed as

where $k_{x}d_{i}=k_{y}d_{i}=2\pi/16$ and $kd_{i}=\sqrt{k_{x}^{2}+k_{y}^{2}}\approx 0.56$, and $\omega=v_{A0}k\left(\tfrac{1}{2}d_{i}k+\sqrt{1+\tfrac{1}{4}d_{i}^{2}k^{2}}\right)$. This intermediate value of $kd_{i}$ is chosen to simultaneously test both the grid-based electron and particle-based ion discretizations. These simulations use $400$ particles/cell initialized with the quasi-quiet start, and apply one pass of conservative binomial smoothing (C). The timestep used is $\Delta t\Omega_{ci}=0.0125$.