Source-linked AI summary
A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations
J. A. Rossmanith, D. C. Seal
TL;DR
The Vlasov system presents computational challenges from high dimensionality and difficulty maintaining conservation and positivity in discretizations. This paper develops a high-order semi-Lagrangian discontinuous Galerkin method with operator splitting, demonstrating accuracy, conservation, positivity, and stability for the 1+1 Vlasov-Poisson system.
Problem
High dimensionality and difficulty exactly maintaining conserved quantities and PDF positivity are central numerical challenges for Vlasov models.
Method
The method combines quasi-1D semi-Lagrangian advection steps, fourth-order operator splitting, high-order discontinuous Galerkin representations, and a modified local DG Poisson solver.
Results
The proposed method is unconditionally stable, fifth-order accurate in space, fourth-order accurate in time, mass conservative, and positivity-preserving.
Takeaways & Limitations
Standard Vlasov-Poisson test cases demonstrated the accuracy and robustness of the proposed scheme.
Takeaways & Limitations
The reported results concern the 1+1 Vlasov-Poisson system; extension to higher-dimensional Vlasov-Poisson equations remains future work.
Abstract
from arXiv · showhide
The Vlasov-Poisson equations describe the evolution of a collisionless plasma, represented through a probability density function (PDF) that self-interacts via an electrostatic force. One of the main difficulties in numerically solving this system is the severe time-step restriction that arises from parts of the PDF associated with moderate-to-large velocities. The dominant approach in the plasma physics community for removing these time-step restrictions is the so-called particle-in-cell (PIC) method, which discretizes the distribution function into a set of macro-particles, while the electric field is represented on a mesh. Several alternatives to this approach exist, including fully Lagrangian, fully Eulerian, and so-called semi-Lagrangian methods. The focus of this work is the semi-Lagrangian approach, which begins with a grid-based Eulerian representation of both the PDF and the electric field, then evolves the PDF via Lagrangian dynamics, and finally projects this evolved field back onto the original Eulerian mesh. In particular, we develop in this work a method that discretizes the 1+1 Vlasov-Poisson system via a high-order discontinuous Galerkin (DG) method in phase space, and an operator split, semi-Lagrangian method in time. Second-order accuracy in time is relatively easy to achieve via Strang operator splitting. With additional work, using higher-order splitting and a higher-order method of characteristics, we also demonstrate how to push this scheme to fourth-order accuracy in time. We show how to resolve all of the Lagrangian dynamics in such a way that mass is exactly conserved, positivity is maintained, and high-order accuracy is achieved. The Poisson equation is solved to high-order via the smallest stencil local discontinuous Galerkin (LDG) approach. We test the proposed scheme on several standard test cases.
1. Introduction
The introduction presents collisionless-plasma modeling challenges and motivates a high-order grid-based semi-Lagrangian DG method for the 1+1 Vlasov-Poisson system.
- Vlasov equations model collisionless plasma dynamics in astrophysical and laboratory settings.
- Numerical challenges: High dimensionality makes the six-dimensional phase-space problem computationally expensive.The system has three spatial and three velocity dimensions, in addition to time.
- Numerical challenges: Exact conservation and positivity of the probability density are difficult to guarantee with many standard Vlasov discretizations.
- Numerical challenges: Moderate-to-large velocities impose severe time-step restrictions relative to the dynamics of interest.
- Existing approaches: Semi-Lagrangian methods evolve a grid-based PDF through Lagrangian dynamics and project it back onto the Eulerian mesh.This retains a grid structure while avoiding small time-step restrictions.
- Current work: The proposed method combines high-order DG representations with operator-split semi-Lagrangian time stepping for the 1+1 Vlasov-Poisson equation.The authors describe it as accurate at relatively low computational expense.
2. Mathematical equations
This section formulates the Vlasov-Poisson model as a kinetic description of a collisionless plasma, specialized to one spatial and one velocity dimension with periodic spatial boundaries.
- The Vlasov system evolves a probability density function in phase space whose moments represent physical quantities such as mass and momentum densities.
- Under non-relativistic and collisionless assumptions, each species PDF obeys an advection equation in position-velocity phase space.
- The paper restricts attention to a single-species Vlasov-Poisson model with electrostatic coupling and a stationary background species.
- Mass, momentum, and additional conserved quantities are used as diagnostics for the numerical discretization.The listed diagnostics include L1 and L2-related quantities and entropy.
- The analyzed problem is the 1+1-dimensional system on [−L, L] × lR with periodic boundary conditions in x.
3. Strang operator splitting
The section explains operator splitting for Vlasov advection: separating spatial and velocity advection enables simple semi-Lagrangian substeps, while higher-order splitting extends temporal accuracy.
- The Vlasov advection equation is split into separate x- and v-advection subproblems with transverse coordinates treated as parameters.
- Each split operator becomes a constant-coefficient advection equation that can be handled with semi-Lagrangian time stepping.The splitting introduces splitting error.
- Strang operator splitting provides a second-order accurate version of the Cheng–Knorr scheme.
- The density and electric field computed in the split algorithm are second-order accurate in time under exact-substep assumptions.
- The present work combines the operator-splitting framework with high-order DG spatial discretization and fourth-order operator splitting.
4. The discontinuous Galerkin (DG) method
The DG formulation represents conservation-law solutions with local polynomial expansions on a Cartesian mesh, using Legendre bases and quadrature to construct high-order approximations.
- The DG method is introduced for a general two-dimensional hyperbolic conservation law on a Cartesian mesh.
- The mesh uses uniform spacing, and the broken finite element space permits elementwise polynomial representations.
- Approximate solutions are expanded in normalized Legendre basis functions on a canonical element.
- For fifth-order spatial accuracy, the representation uses M = 5 and M(M + 1)/2 = 15 Legendre coefficients.
- Initial conditions are obtained by L2 projection, while semidiscrete evolution equations follow from testing the conservation law and integrating by parts.
- The required integrals are evaluated with standard one- and two-dimensional Gaussian quadrature rules.
5. A high-order semi-Lagrangian DG method
This section develops a semi-Lagrangian discontinuous Galerkin method for the 1+1 Vlasov-Poisson system, combining high-order spatial and temporal discretization with mass conservation and positivity preservation. The method advances split advection problems using characteristic shifts and projection, and extends to fourth-order time accuracy.
- The proposed method is unconditionally stable, fifth-order accurate in space, fourth-order accurate in time, mass conservative, and positivity-preserving.
- 5.4. Higher-order splitting: Higher-order operator splitting generalizes the Strang-split method to a fully fourth-order accurate Vlasov-Poisson method.The resulting fourth-order method is summarized as Algorithm 3 after the split solvers and Poisson solver are established.
- 5.1. A toy problem: the 1D advection equation: The basic semi-Lagrangian update exactly advects the initial condition over a time step and projects the result back onto the original mesh.For constant-coefficient advection, the update is unconditionally stable independently of the spatial polynomial order.
- 5.1. A toy problem: the 1D advection equation: The one-dimensional update includes additional ν2 and ν3 terms that ensure stability up to CFL number one.The corresponding stability claim holds for 0 ≤ν ≤1, unlike the related LxW-DG update discussed here.
- 5.2. A semi-Lagrangian DG method for Vlasov-Poisson: The 1+1 Vlasov method applies the one-dimensional advection scheme along Gaussian quadrature lines, shifts each line according to its velocity, and reconstructs the two-dimensional solution.Each line is treated as a one-dimensional polynomial solution, and the reconstructed element solution combines the line solutions with Gaussian quadrature weights.
- 5.2. A semi-Lagrangian DG method for Vlasov-Poisson: The scheme mixes forward and backward characteristic evolution: cell edges move forward to determine projection quadrature points, while solution values are traced backward to the initial time.The forward step is identified as necessary for mass conservation.
5.3. Poisson solver
The Poisson solver uses a local discontinuous Galerkin formulation, rewriting Poisson’s equation as a first-order system and assembling a structured matrix system. Boundary conditions are handled through parameter choices and related formulations.
- The Poisson equation is rewritten as a system of two equations and discretized using the local discontinuous Galerkin method.The formulation follows the LDG approach for solving Poisson equations with discontinuous Galerkin methods.
- One-sided rules are used to evaluate the discrete potential and electric field at element interfaces.
- The resulting matrix equations have a structured form that includes M × M blocks and can be solved efficiently using triangular systems.The matrices A and C can be inverted once at the beginning of the calculation.
- The LDG Poisson formulation can be adapted to Dirichlet boundary conditions by choosing the parameter γ appropriately in the mixed-boundary formulation.Periodic boundary conditions are also handled by carefully choosing β and γ, with β taken as zero without loss of generality.
5.4. High-order split semi-Lagrangian method
The method combines high-order DG spatial discretization with operator-split semi-Lagrangian time integration, extending the approach to fourth-order accuracy while accounting for the electric field’s time dependence.
- Overall method: The complete method is fourth-order accurate in time and fifth-order accurate in space for the Vlasov–Poisson equations.The fourth-order algorithm is summarized in Algorithm 3 and combines the split advection steps with the Poisson solve and electric-field construction.
- Time-dependent electric field: Time-dependent electric fields complicate fourth-order accuracy because Cheng–Knorr splitting does not reduce the system to two constant-coefficient problems.The method constructs temporally accurate electric-field representations rather than treating the field as constant.
- Time-dependent electric field: A fourth-order Taylor polynomial approximates the electric field in time using its value and time derivatives computed from Poisson and moment equations.The second and third time derivatives are expressed through spatial derivatives, which are evaluated using central differences on Legendre coefficients.
- Spatial derivatives: Central finite differences applied directly to DG Legendre coefficients provide the required first- and second-spatial-derivative approximations.The paper verifies these derivative formulas with relative L2 errors for f′(x) and f′′(x) using f(x) = e^sin(2πx).
- Semi-Lagrangian evolution: The resulting scheme uses high-order characteristics for time-dependent advection, transporting interfaces and quadrature points with the time-varying electric field.The semi-Lagrangian DG solver otherwise remains largely unchanged, and the approach requires one Poisson solve per time step.
5.5. Mass conservation and positivity in the mean
The semi-Lagrangian method is proved mass conservative and positive in the mean, with a local limiter extending positivity control to the solution values used by the scheme.
- Mass conservation: Each step of the semi-Lagrangian method conserves total mass.The proof uses periodicity, changes of variables, and index shifting to show that cell contributions preserve the global integral.
- Global positivity: A modified Zhang–Shu limiter enforces non-negativity at test points while remaining local to each element.It leaves already non-negative solutions unchanged and damps high-order corrections when test-point values are negative.
- Positivity in the mean: The positivity theorem assumes that the pre-step solution is non-negative at all points used in the quadrature-based update.This assumption is the key condition in the proof of positivity in the mean.
- Positivity in the mean: Positivity in the mean holds after one coordinate-direction time step independently of the time-step size, provided the prescribed solution-point values are non-negative.The update is represented using shifted solution values and positive quadrature weights.
- Global positivity: The limiter is applied during fourth-order stages and again afterward at tensor-product Gauss–Legendre points to support non-negative integral evaluations.The staged applications preserve positivity in the mean, while the final application assists computations of quantities such as the L1 norm and total energy.
6. Numerical examples
The numerical examples assess the semi-Lagrangian DG scheme on advection, forced Vlasov-Poisson, two-stream instability, and Landau damping problems. They examine convergence, solution structure, positivity, conservation, and agreement with theoretical or previously reported behavior.
- Linear advection: The linear advection test evolves a rotating, compactly supported bump to t = 1 and measures relative L2 errors for Strang and fourth-order splitting.The convergence studies use 5th-order spatial discretization and vary the mesh spacing.
- Forced Vlasov-Poisson: The forced Vlasov-Poisson test compares Strang and fourth-order splitting against an exact solution on multiple phase-space grids.The method uses characteristic evolution with a source-term modification, and the convergence study turns off positivity limiters because they do not guarantee positivity in the mean with a source.
- Two-stream instability: The positivity-preserving limiter removes small negative values that remain in the otherwise reasonable un limited two-stream solution.The same experiment also tracks deviations in the L1 norm, L2 norm, total energy, and entropy.
- Strong Landau damping: Strong Landau damping gives an initial linear decay rate γ1 ≈ −0.292 and a particle-trapping growth rate γ2 = 0.0815.The reported decay rate is closer to Cheng and Knorr's −0.281 than to the other cited comparison values, while the growth rate differs from Heath et al.'s 0.0324.
7. Conclusions and future work
The paper presents a semi-Lagrangian DG method for 1+1 Vlasov-Poisson with high-order accuracy, conservation, positivity, and stability properties. It also identifies higher-dimensional extensions as future work.
- The method is unconditionally stable, fifth-order accurate in space, fourth-order accurate in time, mass conservative, and positivity-preserving.
- The scheme combines quasi-1D semi-Lagrangian advection, fourth-order operator splitting, high-order DG representations, and a modified local DG Poisson solver.
- Standard Vlasov-Poisson test cases demonstrated the proposed scheme’s accuracy and robustness.
- DG represents very rough data effectively, while the semi-Lagrangian framework avoids the decreasing CFL restriction associated with explicit high-order time stepping.
- Future work will extend the approach to higher-dimensional Vlasov-Poisson and to non-relativistic and relativistic Vlasov-Maxwell systems.
Appendix A. Numerical evaluation of conserved quantities
The appendix describes numerical diagnostics for conserved quantities and figure conventions used in strong Landau damping simulations. These diagnostics include norms, energy, and entropy evaluated on specified meshes and resolutions.
- The numerical method uses the L1 norm, L2 norm, total energy, and entropy as diagnostics for conserved quantities.
- Strong Landau damping distribution-function plots use a (128, 256) mesh, fifth-order spatial accuracy, positivity-preserving limiters, and CFL number 2.0.
- Electric-field L2 plots compare (64, 128) and (128, 256) meshes under fifth-order accuracy, positivity preservation, and CFL number 2.0.
- The conserved-quantity evaluations use a functional based on mesh elements, Gauss-Legendre quadrature weights, and quadrature points.
- Strong Landau damping diagnostics plot L1 norm, L2 norm, energy, and entropy across two mesh resolutions.
- At t = 60, both strong and weak Landau damping solutions remain non-zero for |v| > 5, with more activity in the strong-damping case.
Appendix B. Relative L2-norm error in 1D
The appendix defines the relative L2-norm error indicator for one-dimensional DG approximations by comparing an approximate solution with the exact solution over a finite domain.
- The one-dimensional relative L2-norm compares the DG approximation f_h(x) with the exact solution f(x) on x ∈ [a, b].
- The appendix uses an easily computable relative L2-norm indicator based on the total number of grid elements and ∆x = (b−a)/N.
Appendix C. Relative L2-norm error in 2D
The appendix extends the relative L2-norm error formulation to two-dimensional DG approximations over a rectangular domain.
- The two-dimensional relative L2-norm compares f_h(x, y) with f(x, y) on [a_x, b_x] × [a_y, b_y].
- The numerical indicator uses N_x and N_y grid elements with ∆x = (b_x−a_x)/N_x and ∆y = (b_y−a_y)/N_y.