The integration of the vlasov equation in configuration space (2023)

The paper investigates two new use cases for the Boris Spectral Deferred Corrections (Boris-SDC) time integrator for plasma simulations. First, we show that using Boris-SDC as a particle pusher in an electrostatic particle-in-cell (PIC) code can, at least in the linear regime, improve simulation accuracy compared with the standard second order Boris method. In some instances, the higher order of Boris-SDC even allows a much larger time step, leading to modest computational gains. Second, we propose a modification of Boris-SDC for the relativistic regime. Based on an implementation of Boris-SDC in the RUNKO PIC code, we demonstrate for a relativistic Penning trap that Boris-SDC retains its high order of convergence for velocities ranging from $0.5c$ to $>0.99c$.

Simulation of plasmas in electromagnetic fields requires numerical solution of a kinetic equation that describes the time evolution of the particle distribution function. In this paper we propose a finite volume scheme based on integral relation for Poisson brackets to solve the Liouville equation, the most fundamental kinetic equation. The proposed scheme conserves the number of particles, maintains the total-variation-diminishing (TVD) property, and provides high-quality numerical results. Other types of kinetic equations may be also formulated in terms of Poisson brackets and solved with the proposed method including the transport equations describing the acceleration and propagation of Solar Energetic Particles (SEPs), which is of practical importance, since the high energy SEPs produce radiation hazards. The proposed scheme is demonstrated to be accurate and efficient, which makes it applicable to global simulation systems analyzing space weather.

We study the conservation properties of the Hermite-discontinuous Galerkin (Hermite-DG) approximation of the Vlasov-Maxwell equations. In this semi-discrete formulation, the total mass is preserved independently for every plasma species. Further, an energy invariant exists if central numerical fluxes are used in the DG approximation of Maxwell's equations, while a dissipative term is present when upwind fluxes are employed. In general, traditional temporal integrators might fail to preserve invariants associated with conservation laws during the time evolution. Hence, we analyze the capability of explicit and implicit Runge-Kutta (RK) temporal integrators to preserve such invariants. Since explicit RK methods can only ensure preservation of linear invariants but do not provide any control on the system energy, we consider modified explicit RK methods in the family of relaxation Runge-Kutta methods (RRK). These methods can be tuned to preserve the energy invariant at the continuous or semi-discrete level, a distinction that is important when upwind fluxes are used in the discretization of Maxwell's equations since upwind provides a numerical source of energy dissipation that is not present when central fluxes are used. We prove that the proposed methods are able to preserve the energy invariant and to maintain the semi-discrete energy dissipation (if present) according to the discretization of Maxwell's equations. An extensive set of numerical experiments corroborates the theoretical findings. It also suggests that maintaining the semi-discrete energy dissipation when upwind fluxes are used leads to an overall better accuracy of the method relative to using upwind fluxes while forcing exact energy conservation.

Direct kinetic solvers enable high-fidelity simulation of multiscale plasmas in a number of applications including space physics and propulsion, materials processing, astrophysics, and nuclear fusion. While they eliminate the statistical noise associated with particle-in-cell methods, they are associated with higher dimensionalities. Thus, a detailed understanding of grid-point requirements is required to design efficient meshes so that direct kinetic solvers can be feasibly employed in general settings without compromising on predictivity. The grid-point requirements of a direct kinetic solver employing the Vlasov-Poisson-BGK equations are characterized using an electrostatic and weakly collisional plasma plume expansion model problem, which is unsteady and spatially inhomogeneous with significant deviations from equilibrium. It is demonstrated that at least two to four points per the appropriate Debye length and thermal velocity are necessary to resolve macroscopic density gradients and thus the lowest-order macroscopic quantities, with more stringent requirements for higher-order quantities. Local charge separation and the distribution function itself require at least an additional order of magnitude in resolution for comparable accuracy, as they require the resolution of gradients in the distribution function. Collisions impede plume expansion and introduce secondary flow and field structures, but do not significantly relax the grid requirements despite their smoothing action in velocity space. While specific numerical requirements necessarily depend on the solver, plasma configuration, and collision model, trends in the variation of the elucidated grid-point requirements with the quantity of interest and plasma collisionality can be generalizable across problems with comparable physics. Knowledge of similarly derived trends can contribute to efficient direct kinetic simulation of unsteady and spatially inhomogeneous plasmas of a variety of configurations and collisionalities.

In this study, we design a specific geometry of plasma particles and further translate it into mathematical form, i.e., formulated time-fractional semi relativistic Vlasov Maxwell system. This extended version of model is very significant and ground-breaking because it has ability to study the behavior of plasma particles at macroscopic and microscopic time-evaluation scales, which is not studied yet. This model is further tackled with numerical strategy, projected in accordance with spectral and finite difference approximations. The numerical results demonstrate that there are certain variations of the plasma particles that were out of sight at fractional scale and dynamically, we have entered into the real profundity of the problem. Numerical convergence of the projected technique is also inspected. These designed tools, such as the plasma and numerical methods, are very significant and can say that it is a remarkable contribution in this field. In addition, the established numerical structure can utilize to scrutinize the required numerical solution of highly non-linear and fractional problems.

A Fourier transformation based unified gas-kinetic scheme (UGKS) is proposed to simulate the kinetic behaviors of plasma in cylindrical coordinates $(r,\theta ,v_r,v_{\theta })$, the model is depicted by Vlasov–Poisson equations coupled with Bhatnagar–Gross–Krook(BGK) collision term. Fourier transformation is applied to the equations for $\theta$ and a series of equations are gotten. Based on Strang-splitting strategy, Vlasov equations are divided into two parts, the transport-collision part solved by a multiscale gas-kinetic scheme, and acceleration part solved by Runge–Kutta method. The algorithm is applied on charge separation problem at plasma edge and Z-pinch configuration. Numerical results show our scheme can give a clear kinetic picture of plasma in 2D phase space, and also can capture the process from non-equilibrium to equilibrium state by Coulomb collisions.

Journal of Computational Physics, Volume 267, 2014, pp. 7-27

Semi-Lagrangian schemes with various splitting methods, and with different reconstruction/interpolation strategies have been applied to kinetic simulations. For example, the order of spatial accuracy of the algorithms proposed in Qiu and Christlieb (2010) [29] is very high (as high as ninth order). However, the temporal error is dominated by the operator splitting error, which is second order for Strang splitting. It is therefore important to overcome such low order splitting error, in order to have numerical algorithms that achieve higher orders of accuracy in both space and time. In this paper, we propose to use the integral deferred correction (IDC) method to reduce the splitting error. Specifically, the temporal order accuracy is increased by r with each correction loop in the IDC framework, where $r=1,2$ for coupling the first order splitting and the Strang splitting, respectively. The proposed algorithm is applied to the Vlasov–Poisson system, the guiding center model, and two dimensional incompressible flow simulations in the vorticity stream-function formulation. We show numerically that the IDC procedure can automatically increase the order of accuracy in time. We also investigate numerical stability of the proposed algorithm via performing Fourier analysis to a linear model problem.

Computer Physics Communications, Volume 254, 2020, Article 107351

In this paper, our goal is to efficiently solve the Vlasov equation on GPUs. A semi-Lagrangian discontinuous Galerkin scheme is used for the discretization. Such kinetic computations are extremely expensive due to the high-dimensional phase space. The SLDG code, which is publicly available under the MIT license, abstracts the number of dimensions and uses a shared codebase for both GPU and CPU based simulations. We investigate the performance of the implementation on a range of both Tesla (V100, Titan V, K80) and consumer (GTX 1080 Ti) GPUs. Our implementation is typically able to achieve a performance of approximately 470 GB/s on a single GPU and 1600 GB/s on four V100 GPUs connected via NVLink. This results in a speedup of about a factor of ten (comparing a single GPU with a dual socket Intel Xeon Gold node) and approximately a factor of 35 (comparing a single node with and without GPUs). In addition, we investigate the effect of single precision computation on the performance of the SLDG code and demonstrate that a template based dimension independent implementation can achieve good performance regardless of the dimensionality of the problem.

Journal of Computational Physics, Volume 288, 2015, pp. 66-85

Computer Physics Communications, Volume 191, 2015, pp. 65-73

Particle-in-cell merging algorithms aim to resample dynamically the six-dimensional phase space occupied by particles without distorting substantially the physical description of the system. Whereas various approaches have been proposed in previous works, none of them seemed to be able to conserve fully charge, momentum, energy and their associated distributions. We describe here an alternative algorithm based on the coalescence of N massive or massless particles, considered to be close enough in phase space, into two new macro-particles. The local conservation of charge, momentum and energy are ensured by the resolution of a system of scalar equations. Various simulation comparisons have been carried out with and without the merging algorithm, from classical plasma physics problems to extreme scenarios where quantum electrodynamics is taken into account, showing in addition to the conservation of local quantities, the good reproducibility of the particle distributions. In case where the number of particles ought to increase exponentially in the simulation box, the dynamical merging permits a considerable speedup, and significant memory savings that otherwise would make the simulations impossible to perform.

Computer Physics Communications, Volume 198, 2016, pp. 47-58

We describe a spectral method for the numerical solution of the Vlasov–Poisson system where the velocity space is decomposed by means of an Hermite basis, and the configuration space is discretized via a Fourier decomposition. The novelty of our approach is an implicit time discretization that allows exact conservation of charge, momentum and energy. The computational efficiency and the cost-effectiveness of this method are compared to the fully-implicit PIC method recently introduced by Markidis and Lapenta (2011) and Chen etal. (2011). The following examples are discussed: Langmuir wave, Landau damping, ion-acoustic wave, two-stream instability. The Fourier–Hermite spectral method can achieve solutions that are several orders of magnitude more accurate at a fraction of the cost with respect to PIC.

Journal of Computational Physics, Volume 340, 2017, pp. 470-497

We have developed a deterministic conservative solver for the inhomogeneous Fokker–Planck–Landau equation coupled with the Poisson equation, which is a classical mean-field primary model for collisional plasmas. Two subproblems, i.e. the Vlasov–Poisson problem and homogeneous Landau problem, are obtained through time-splitting methods, and treated separately by the Runge–Kutta Discontinuous Galerkin method and a conservative spectral method, respectively. To ensure conservation when projecting between the two different computing grids, a special conservation routine is designed to link the solutions of these two subproblems. This conservation routine accurately enforces conservation of moments in Fourier space. The entire numerical scheme is implemented with parallelization with hybrid MPI and OpenMP. Numerical experiments are provided to study linear and nonlinear Landau Damping problems and two-stream flow problem as well.