Source-linked AI summary
A low-rank projector-splitting integrator for the Vlasov--Poisson equation
Lukas Einkemmer, Christian Lubich
TL;DR
Direct grid discretization of Vlasov–Poisson is costly in high-dimensional phase space. The paper applies dynamical low-rank approximation with projector splitting, including a hierarchical variant, and reports that small ranks capture the physics of several simulations while substantially reducing computational cost.
Problem
Direct discretization of Vlasov–Poisson in high-dimensional phase space requires O(n^2d) storage and can be computationally prohibitive.
Source: §1 · p. 1Method
The method constrains Vlasov–Poisson dynamics to a low-rank manifold, splits the tangent-space projection, and can use a hierarchical low-rank approximation.
Source: §1 · p. 3Results
Small ranks capture the physics of linear Landau damping, a two-stream instability, and a plasma echo, substantially reducing computational cost.
Source: §1 · p. 3Takeaways & Limitations
The scheme stays within the low-rank manifold without forming and truncating high-rank tensors, helping reduce storage and computational cost.
Source: §1 · p. 3With the electric field assumed fixed during the relevant step, naive Strang splitting still gives only a first-order method.
Source: §2.3 · p. 9
Abstract
from arXiv · showhide
Many problems encountered in plasma physics require a description by kinetic equations, which are posed in an up to six-dimensional phase space. A direct discretization of this phase space, often called the Eulerian approach, has many advantages but is extremely expensive from a computational point of view. In the present paper we propose a dynamical low-rank approximation to the Vlasov--Poisson equation, with time integration by a particular splitting method. This approximation is derived by constraining the dynamics to a manifold of low-rank functions via a tangent space projection and by splitting this projection into the subprojections from which it is built. This reduces a time step for the six- (or four-) dimensional Vlasov--Poisson equation to solving two systems of three- (or two-) dimensional advection equations over the time step, once in the position variables and once in the velocity variables, where the size of each system of advection equations is equal to the chosen rank. By a hierarchical dynamical low-rank approximation, a time step for the Vlasov--Poisson equation can be further reduced to a set of six (or four) systems of one-dimensional advection equations, where the size of each system of advection equations is still equal to the rank. The resulting systems of advection equations can then be solved by standard techniques such as semi-Lagrangian or spectral methods. Numerical simulations in two and four dimensions for linear Landau damping, for a two-stream instability and for a plasma echo problem highlight the favorable behavior of this numerical method and show that the proposed algorithm is able to drastically reduce the required computational effort.
1. Introduction.
The paper adapts dynamical low-rank approximation and projector splitting to reduce the cost of Vlasov–Poisson simulations in high-dimensional phase space. It reports that small ranks capture the physics in several test problems while substantially reducing computational cost.
1. Introduction.
- Direct Eulerian discretization of the up-to-six-dimensional phase space requires O(n^2d) storage and can be prohibitive.Here n is the number of grid points per direction.Source: §1 · p. 1
- Filamentation creates progressively smaller phase-space scales, while the linear part of the equation can nevertheless be represented exactly with rank two.The paper also notes that kinetic equilibria formed by linear combinations of Gaussians can be represented efficiently at low rank.Source: §1 · pp. 2–3
- Sparse-grid approaches face large mixed derivatives and can struggle with Gaussians, and their cited use can develop instabilities when too few degrees of freedom are used.A cited four-dimensional Landau-damping scheme reduced memory by approximately a factor of 10 and runtime by a factor of 2, but remedied its Gaussian deficiency with a δf scheme that performs poorly when particle density deviates significantly from that form.Source: §1 · p. 2
- The method constrains Vlasov–Poisson dynamics to a low-rank manifold and splits the tangent-space projection into a sequence of lower-dimensional advection problems.The first-order scheme uses rank-r systems in x and v, with an intermediate system of r^2 ordinary differential equations.Source: §1 · p. 3
- A hierarchical Tucker approximation can further reduce a time step to one-dimensional advection equations and integrals.The paper adapts a projector-splitting integrator for low-rank Tucker tensors.Source: §1 · p. 3
- The scheme stays within the low-rank manifold without forming and truncating a high-rank tensor, helping reduce storage and computational cost.Depending on the problem, physical or velocity space can be solved directly while retaining the hierarchical approximation elsewhere.Source: §1 · p. 3
- Small-rank simulations capture the physics of linear Landau damping, a two-stream instability, and a plasma echo while substantially reducing computational cost.The introduction presents these simulations as evidence for the proposed method.Source: §1 · p. 3
2. Description of the numerical algorithm.
The algorithm constrains the Vlasov–Poisson dynamics to a low-rank manifold and splits its tangent-space projection into sequential updates. Spectral and semi-Lagrangian methods solve the resulting equations, while a midpoint electric-field estimate makes a second-order scheme possible.
2. Description of the numerical algorithm.
- The approximation represents f using r time-dependent x- and v-basis functions coupled by a coefficient matrix, with orthonormal bases.Here r is the rank of the representation.Source: §2 · p. 4
- Splitting the tangent-space projector into three terms yields sequential evolution equations for the low-rank approximation.The first-order method uses Lie–Trotter splitting; a second-order extension uses Strang splitting.Source: §2 · p. 5
2.1. Spectral method.
- Fourier transformation makes the linear part the only source of stiffness, so an exponential integrator can integrate it exactly and remove the CFL condition.For the velocity-space update, the electric field is frozen at the start of the step, retaining a first-order approximation.Source: §2.1 · p. 8
2.2. Semi-Lagrangian method.
- Diagonalizing the symmetric coefficient matrix turns the semi-Lagrangian update into constant-coefficient one-dimensional advection equations.After solving them, the coordinate transformation is undone; the procedure is repeated across spatial directions.Source: §2.2 · p. 9
2.3. Second-order integrator.
- A naive Strang splitting remains first order because the electric field is fixed during one update; an embedded half-step estimates the midpoint field for the second-order scheme.The method uses that estimate in the middle Strang step and composes two Lie steps in opposite order.Source: §2.3 · p. 9
3. A hierarchical low-rank algorithm.
The hierarchical algorithm adds low-rank structure within the spatial and velocity factors to reduce the dimensions of the advection equations further. The paper describes the construction for four dimensions and extends it to six dimensions.
3. A hierarchical low-rank algorithm.
- The hierarchical approach further reduces the d-dimensional problems from the preceding algorithm when the spatial and velocity factors also admit low-rank structure.The paper first develops the four-dimensional case and then states the six-dimensional extension.Source: §3 · p. 9
3.1. The hierarchical low-rank algorithm in four dimensions.
- In four dimensions, separate spatial and velocity ranks allow the representation to capture cases where only velocity-space structure is low rank.The resulting one-variable evolution equations can use spectral or semi-Lagrangian methods.Source: §3.1 · pp. 10, 15
- In four dimensions, evolution equations cost O(nR2) operations and O(nR) storage; coefficient computation costs O(nR2 + R4) operations.R = max(r, rx, rv, rE).Source: §3.1 · p. 15
3.2. The hierarchical low-rank algorithm in six dimensions.
- The six-dimensional extension updates three spatial coordinate factors and then the core, using one-dimensional evolution equations.The section states that its derivation is similar to the four-dimensional case and gives the relevant evolution equations.Source: §3.1 · p. 16, §3.2 · p. 16
4. Numerical results.
The numerical tests examine low-rank and hierarchical low-rank splitting methods on Landau damping, two-stream instability, and plasma echo problems. Across these cases, low ranks reproduce key dynamics, while the echo simulations show substantial computational-effort reduction compared with direct simulation.
4.1. Linear Landau damping.
- Rank r = 5 closely matches the analytic Landau-damping solution, with energy error around 10^-8 and mass and L2 errors indistinguishable from machine precision.Source: §4.1 · p. 19, Figure 4.1 · p. 20
- In four-dimensional aligned Landau damping, rank (5, 5, 5) predicts the decay rate, while rank (10, 10, 10) reduces electric energy to 10^-8, about two orders better.Mass, energy, and L2 errors remain below 10^-6 in all configurations.Source: Figure 4.2 · p. 20, §4.1 · p. 21
- For non-aligned four-dimensional Landau damping, ranks (10, 10, 10) and (15, 15, 15) reduce electric energy to 5 · 10^-6 and 5 · 10^-7, respectively.The rank r can be smaller than rx and rv while still producing good results; mass, energy, and L2 errors stay below 10^-5.Source: Figure 4.3 · p. 21, §4.1 · p. 21
4.2. Two-stream instability.
- In the two-dimensional two-stream instability, rank r = 5 agrees excellently with direct Eulerian results in the linear regime, while ranks from 10 keep nonlinear-regime electric energy approximately constant.Most mass and energy error accrues during exponential electric-energy growth, with no observable long-term invariant drift; L2 is conserved to machine precision.Source: Figure 4.4 · p. 22, §4.2 · p. 22
- For four-dimensional two-stream instability, rank (5, 5, 5) gives excellent linear-regime agreement, and rank (10, 10, 10) improves mass, energy, and L2 errors by about one order of magnitude.Source: Figure 4.5 · p. 23, §4.2 · p. 23
4.3. Plasma echo.
- At r = 5, both plasma echoes are resolved; r = 10 resolves secondary-echo growth and decay, with invariant errors below 10^-10 and substantially reduced effort versus direct simulation.The simulation uses 4096 velocity-direction and 512 space-direction grid points.Source: Figure 4.6 · p. 24, §4.3 · p. 24