Discontinuous Galerkin algorithms for fully kinetic plasmas
J. Juno, A. Hakim, J. TenBarge, E. Shi, and W. Dorland
TL;DR
This paper introduces a high-order discontinuous Galerkin algorithm for solving the Vlasov-Maxwell system in plasma physics, emphasizing efficiency and accuracy in high-dimensional simulations.
Contribution
It develops a novel high-order discontinuous Galerkin method with cost-reduction features for kinetic plasma simulations, validated through various benchmarks.
Findings
Accurate high-order solutions for plasma distribution functions.
Effective reduction of computational cost in high-dimensional simulations.
Successful benchmarks demonstrating method efficacy.
Abstract
We present a new algorithm for the discretization of the Vlasov-Maxwell system of equations for the study of plasmas in the kinetic regime. Using the discontinuous Galerkin finite element method for the spatial discretization, we obtain a high order accurate solution for the plasma's distribution function. Time stepping for the distribution function is done explicitly with a third order strong-stability preserving Runge-Kutta method. Since the Vlasov equation in the Vlasov-Maxwell system is a high dimensional transport equation, up to six dimensions plus time, we take special care to note various features we have implemented to reduce the cost while maintaining the integrity of the solution, including the use of a reduced high-order basis set. A series of benchmarks, from simple wave and shock calculations, to a five dimensional turbulence simulation, are presented to verify the…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40| Lagrange Tensor | Polynomial Order | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Dimension | ||||||||
| 2 | 4 | 9 | 16 | 25 | 36 | 49 | 64 | |
| 3 | 8 | 27 | 64 | 125 | 216 | 343 | 512 | |
| 4 | 16 | 81 | 256 | 625 | 1296 | 2401 | 4096 | |
| 5 | 32 | 243 | 1024 | 3125 | 7776 | 16807 | 32768 | |
| 6 | 64 | 729 | 4096 | 15625 | 46656 | 117649 | 262144 |
| Serendipity Element | Polynomial Order | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Dimension | ||||||||
| 2 | 4 | 8 | 12 | 17 | 23 | 30 | 38 | |
| 3 | 8 | 20 | 32 | 50 | 74 | 105 | 144 | |
| 4 | 16 | 48 | 80 | 136 | 216 | 328 | 480 | |
| 5 | 32 | 112 | 192 | 352 | 592 | 952 | 1472 | |
| 6 | 64 | 256 | 448 | 880 | 1552 | 2624 | 4256 |
| Base Polynomial Order | 1 | 2 | 3 | 4 | 5 | |
| Dimension | ||||||
| 1X3V | 16 | 108 | 320 | 875 | 1728 | |
| 2X3V | 32 | 432 | 1600 | 6125 | 13824 | |
| 3X3V | 64 | 1728 | 8000 | 42875 | 110592 |
| Cost(Original/New) | Base Polynomial Order | 1 | 2 | 3 | 4 | 5 |
| Dimension | ||||||
| 1X3V | 1 | |||||
| 2X3V | 1 | |||||
| 3X3V | 1 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Discontinuous Galerkin algorithms for fully kinetic plasmas
J. Juno
A. Hakim
J. TenBarge
E. Shi
W. Dorland
IREAP, University of Maryland, College Park, MD, 20742
Princeton Plasma Physics Laboratory, Princeton, NJ 08543-0451
Princeton University, Princeton, NJ, 08544
Abstract
We present a new algorithm for the discretization of the non-relativistic Vlasov-Maxwell system of equations for the study of plasmas in the kinetic regime. Using the discontinuous Galerkin finite element method for the spatial discretization, we obtain a high order accurate solution for the plasma’s distribution function. Time stepping for the distribution function is done explicitly with a third order strong-stability preserving Runge-Kutta method. Since the Vlasov equation in the Vlasov-Maxwell system is a high dimensional transport equation, up to six dimensions plus time, we take special care to note various features we have implemented to reduce the cost while maintaining the integrity of the solution, including the use of a reduced high-order basis set. A series of benchmarks, from simple wave and shock calculations, to a five dimensional turbulence simulation, are presented to verify the efficacy of our set of numerical methods, as well as demonstrate the power of the implemented features.
keywords:
Discontinuous Galerkin , Vlasov-Maxwell
††journal: Journal of Computational Physics
Contents
1 Introduction
Plasmas are ubiquitous in nature, and the study of plasmas has application to a wide variety of problems, from the development of nuclear fusion, to understanding the dynamic interaction between the solar wind and the Earth’s magnetosphere, to elucidating the mysteries of large scale astrophysical phenomena such as binary star collisions or the accretion disks of black holes. Unfortunately, many plasmas of interest are only weakly collisional and far from equilibrium, making the system best described by kinetic theory, in this case, the Vlasov equation. The use of kinetic theory significantly complicates the theoretical analysis and simulation of the plasma’s dynamics due to the increased dimensionality of the corresponding equations, which are solved in a combined position and velocity phase space, along with the large collection of waves and instabilities that the kinetic system supports. This complication is usually mitigated in one of two ways: either a reduction of the system via asymptotic expansions of the corresponding equations in appropriate limits for the problem of interest or direct numerical simulation of the plasma by approximating it as a collection of “macroparticles,” and employing the particle-in-cell (PIC) algorithm [1]. At scales much larger than the ion gyroradius, fluid models for the plasma are excellent tools for analyzing the macroscopic evolution of the plasma. Likewise, reductions of the system which retain kinetic effects, such as the gyrokinetic expansion of the Vlasov equation [2, 3], which averages over the fast cyclotron motion of the particles in the plasma to eliminate one velocity space dimension, have been quite lucrative in both fusion and astrophysical contexts. The PIC algorithm has a long history of success in revealing many novel kinetic features of the plasma’s dynamics because, outside of the approximation of integrating along “macroparticle” trajectories as opposed to individual particle trajectories, the algorithm makes no approximations of the plasma’s dynamics and retains the vast array of kinetic phenomena intrinsic to the system.
However, these two approaches are not without their disadvantages. Fluid models derived from asymptotic expansions of the Vlasov equation often require sub-grid models for the microphysics being approximated, and, while physically motivated, these closures are not always well-suited for the plasma being modeled. Even kinetic reductions which are closed, sometimes rigorously so, such as the aforementioned gyrokinetics, are almost certainly not valid in all contexts; for example, large amplitude fluctuations and waves above the cyclotron frequency observed in the solar wind and solar corona are beyond the scope of traditional gyrokinetic formulations. Thus, a tool like the PIC algorithm is ideal for studying physics in regimes where no tractable and physically reasonable asymptotic reduction exists, but the approximation of the plasma as a collection of macroparticles introduces counting noise into the resulting solutions, which can pollute the results quite severely, or at the very least complicate the analysis. One can always mollify this concern by increasing the number of particles in the simulation, but the counting noise decreases like , where is the number of particles per grid cell. This scaling can make the study of physics where the signal-to-noise ratio is low especially challenging. For example, Camporeale et al. [4] have demonstrated that a large number of particles-per-cell is required to correctly identify wave-particle resonances and compare well with linear theory. One could employ the delta-f PIC method [5], but noise mitigation techniques such as the delta-f PIC method can break down if the distribution function deviates significantly from its initial value.
Therefore, if the problem of interest requires a full Vlasov-Maxwell description and is characterized by a low signal-to-noise ratio, an alternative approach which directly solves the Vlasov equation and gives access to a noise-free solution is desired. While computationally infeasible in the past due to the need to solve a six-dimensional system (three position and three velocity), plus time, to accurately capture the dynamics of the plasma, in recent years, directly solving the Vlasov equation has become a more popular alternative approach to illuminating the microphysics of the plasma’s evolution. Previous Vlasov work has focused almost exclusively on the so-called hybrid framework, treating the electrons as a fluid, usually massless and isothermal, to save significantly on computational cost [6, 7, 8, 9, 10, 11, 12, 13, 14], with some exceptions [15, 16]. In this paper, we present the addition of a multi-species Vlasov-Maxwell solver to Gkeyll, a modular framework in which a variety of plasma physics and fluid dynamics solvers are currently being built [17, 18, 19, 20, 21, 22, 23, 24]. Within Gkeyll, the Vlasov-Maxwell system is discretized in space with a discontinuous Galerkin finite element method, and discretized in time with a strong stability preserving Runge-Kutta method to form a fully explicit update. The discontinuous Galerkin method combines the power of finite element methods, including high order accuracy and the ability to handle complicated geometries, with the advantages of finite volume methods, including the introduction of limiters to enforce stability and physicality of the solution, and locality of data, for efficient parallelization [25, 26, 27]. The high order accuracy and locality of data in particular make our approach especially advantageous; higher order polynomials provide a level of accuracy equivalent to refining the grid at a fraction of the cost, and the locality of data significantly reduces the amount of communication required in the update, enabling the algorithm to scale to well on both standard computing architecture and many-core devices. In fact, the discontinuous Galerkin algorithm has been gaining increased attention as a means of discretizing high dimensional transport equations and various flavors of the Vlasov equation, including Vlasov-Poisson, Vlasov-Ampere, and the aforementioned Vlasov-Maxwell [28, 29, 30].
The paper is organized as follows: Section 2 provides an overview of the relevant plasma kinetic equation, with discussion of relevant conservation properties in the continuous system. Section 3 describes the Runge-Kutta discontinuous Galerkin discretization and proves relevant conservation properties of the discrete system in the continuous time limit. Section 4 covers the details of the implementation of the numerical algorithm, including the choice of basis functions to mitigate the computational cost in higher dimensions. Section 5 demonstrates the functionality of the algorithm, including scaling studies and performance analysis of the code on high dimensional benchmarks. Finally, Section 6 summarizes our findings and provides future directions for further numerical improvements while discussing current problems in plasma physics which are now within grasp with this new numerical tool.
2 The Vlasov-Maxwell system
2.1 Basic equations
The time evolution for the distribution function for species in a plasma is given by the Vlasov equation,
[TABLE]
where is the Lorentz force. In the Lorentz force, and are the charge and mass of species respectively, and and are the electric and magnetic fields. The operators and are the gradient operators in configuration and velocity space respectively. Here we are considering two simultaneous limits: the non-relativistic111The special relativistic limit of the Vlasov equation can be obtained by substitution of the Lorentz boost factors where appropriate and a change of variables from velocity to momentum,
(2)
where is the Lorentz boost factor,
(3)
and the collisionless222The weakly collisional or collisional limits can be obtained by the addition of the Landau-Fokker-Planck [31] collision operator to the right hand side of Eq. (1),
(4)
where
(5)
is the collision frequency of species colliding with species , and is the Landau tensor,
(6)
Here, is the identity tensor. In the definition of the collision frequency, is the Debye length. Often, the term in the logarithm, , is abbreviated as and called the plasma parameter. Reductions of the collision operator, such as the Bhatnagar-Gross-Krook (BGK) [32] collision operator,
(7)
where is some equilibrium distribution function to which the distribution function should relax, are also sometimes considered the weakly collisional or collisional limit. However, given the form of the Landau-Fokker-Planck collision operator, reduced collision operators which maintain a Fokker-Planck structure, i.e., a drag and diffusion term, such as the Lenard-Bernstein [33, 34, 35] collision operator,
(8)
where and are related to the first and second moments of the distribution function respectively, are often preferred. Note that, while the Landau-Fokker-Planck operator naturally generalizes to a multi-species plasma, as well as the relativistic Vlasov equation, reduced collision operators such as the BGK and Lenard-Bernstein collision operators may require significant modification to include the effects of cross-species collisions and relativistic effects. limits.
Since the motion of charged particles creates currents and electromagnetic fields, the electric and magnetic fields in the Vlasov equation evolve self-consistently by the coupling of the Vlasov equation to Maxwell’s equations,
[TABLE]
where and are the charge and current densities of the plasma. Charge and current densities are determined by moments of the distribution function
[TABLE]
and
[TABLE]
where the moment operator for any function is defined as
[TABLE]
If there are external fields, and , e.g., created by coils or electrodes, we can replace and in the Lorentz force term. External fields do not appear in Maxwell’s equations.
2.2 Conservation Properties of the Vlasov-Maxwell system
As we mentioned previously, we are interested in deriving a discretization of the non-relativistic, collisionless, Vlasov-Maxwell system. Henceforth, when referring to the Vlasov equation, or the Vlasov-Maxwell system of equations, we emphasize that we are referring to the form of the Vlasov equation given in in Eq. (1). Before proceeding to the derivation of the discretization, there are several key conservation properties of the continuous Vlasov-Maxwell system which we review here for the purposes of exploring these same conservation relations in the discrete system. To prove various properties of the Vlasov-Maxwell system we need to assume that faster than for finite . Also, we will assume that either the configuration space is periodic, or that the distribution function vanishes on the boundaries.
Proposition 1**.**
The Vlasov-Maxwell system conserves particles.
Proof.
Upon integration of Eq. (1) over all of phase-space and summation over all species, we obtain
[TABLE]
This relation is a straightforward consequence of the fact that the Vlasov-Maxwell system is written as a conservation law in phase-space. We further note that this relation holds independently for each species because in the absence of recombination and ionization, there is no mechanism for conversion of one particle species to another. ∎
Proposition 2**.**
The collisionless Vlasov-Maxwell system conserves the norm of the distribution function, i.e.,
[TABLE]
where the integration is taken over the complete phase-space, .
Proof.
The Vlasov-Maxwell system can be written as a non-linear advection equation in phase-space. To do this, we introduce the phase-space velocity vector and the phase-space gradient operator . The flow in phase-space is incompressible, i.e, , because is the velocity coordinate and thus has no configuration space dependence, and likewise, the Lorentz force has no velocity space divergence. The latter point is slightly subtle because while the electric field has no velocity space dependence, the force has no velocity space divergence by properties of the cross product. Eq. (1) can then be written as
[TABLE]
where is the distribution function. If we multiply Eq. (18) by and integrate over all of phase-space we obtain
[TABLE]
We have integrated by parts and used boundary conditions to eliminate the surface term. Now, using incompressibility of phase-space, we can write . Using this relation in the above expression and converting the volume integral into a surface integral, combined with the use of boundary conditions, gives the desired conservation law. ∎
Proposition 3**.**
The collisionless Vlasov-Maxwell system conserves the entropy of the system333The inclusion of the minus sign in the definition of the entropy is in concordance with the traditional physics definition of the entropy. It is common in the theory of hyperbolic conservation laws to drop the minus sign. This definition has the effect that the physicist’s entropy is a non-decreasing quantity, while the mathematician’s entropy is a non-increasing quantity.,
[TABLE]
Proof.
We multiply Eq. (18) by , integrate over all of phase space, and after a bit of algebra, we obtain,
[TABLE]
In the first two terms, we have used the fact that by the product rule, and in the last term, we have used integration by parts and eliminated the surface term. Now, since , we can see that the latter two terms just become the Vlasov equation written in the form of Eq. (18). Since in the collisionless limit the Vlasov equation equals 0, the latter two terms vanish and we are left with our desired conservation relation. ∎
Proposition 4**.**
The Vlasov-Maxwell system conserves total (particles plus field) momentum.
Proof.
We multiply Eq. (1) by , integrate over all of phase-space, and sum over all species to obtain,
[TABLE]
In the last term, we have used integration by parts and . To make further progress, we consider Maxwell’s equations. Taking the cross-product of Eq. (9) with , and the cross-product of Eq. (10) with and subtracting the resulting equations we obtain
[TABLE]
Now, for any vector field we have and . Using these vector identities and the divergence Eqns. (11) and (12) gives
[TABLE]
Finally, inserting Eq. (24) into Eq. (22), integrating by parts, and using configuration space boundary conditions gives
[TABLE]
The first term is the particle momentum, and the second term is the field momentum. ∎
Proposition 5**.**
The Vlasov-Maxwell system conserves total (particles plus field) energy.
Proof.
We multiply Eq. (1) by , integrate over all phase-space, and sum over all species to obtain,
[TABLE]
where we have used integration by parts and to evaluate the third term on the left hand side. To make further progress, we again examine Maxwell’s equations. Taking the dot product of Eq. (10) with and the dot product of Eq. (9) with and adding the resulting equations gives us
[TABLE]
Using this result in Eq. (26), integrating by parts, and using configuration space boundary conditions gives the total energy conservation law
[TABLE]
The first term in this expression is the total particle energy, and the second two terms make up the electromagnetic field energy. ∎
3 The Discrete Vlasov-Maxwell system
3.1 The semi-discrete Vlasov equation
To derive the semi-discrete Vlasov equation using a discontinuous Galerkin algorithm, we begin with the Vlasov equation in the form of Eq. (18). To discretize this equation, we introduce a phase-space mesh with cells , and introduce the following piecewise polynomial approximation space for the distribution function
[TABLE]
where is some space of polynomials. The problem can now be stated as finding such that, for all ,
[TABLE]
for all test functions . In this discrete weak-form, is a numerical flux function, where is an outward unit vector on the surface of the cell . Further, the subscript indicates the discrete solution and the notation () indicates that the function is evaluated just inside (outside) the location on the surface . The particular polynomial space we choose is the so-called Serendipity element space [36]. With the choice of polynomial space, we can then say the discrete solution is,
[TABLE]
where are a set of polynomials chosen such that they lie in the aforementioned Serendipity element space. For the numerical flux, we use the local Lax-Friedrichs flux,
[TABLE]
where . This particular numerical flux function is a common choice for scalar hyperbolic equations due to its simplicity and efficiency. An important note is that the Lax-Friedrichs flux can be reduced to a pure central flux function by setting , where the average of the cell interface values is used to evaluate the surface integrals. A similar reduction can be done with the appropriate choice of to create the standard upwind numerical flux, which uses the sign of to determine the direction of the flow of the distribution function across a cell interface. The surface and the volume integrals in Eq. (30) are replaced by Gaussian quadrature of an appropriate order to ensure that the discrete integrals are performed exactly.
3.2 The semi-discrete Maxwell system
For Maxwell’s equations, we also discretize the system with a discontinuous Galerkin algorithm. We denote the restriction of the phase-space mesh, , to configuration space by . The cells in configuration space are denoted by , for , where are the number of configuration space cells. We introduce the solution space
[TABLE]
i.e, the restriction of the set to configuration space .
While deriving the discrete weak-form, we need to evaluate volume integrals which include terms of the form , for . For example, consider
[TABLE]
Gauss’ law can then be used to convert one volume integral into a surface integral
[TABLE]
where is the (vector) area-element that points in the direction of the outward normal to the configuration space cell . Using these expressions, we can now write the discrete weak-form of Maxwell’s equations as
[TABLE]
and
[TABLE]
for all . Here and are, respectively, the values of the electric and magnetic field at the cell interfaces.
We consider two methods of obtaining the cell interface fields needed in the discrete weak-form of Maxwell’s equations: central fluxes and upwind fluxes. As we will see later, both numerical flux functions have advantages and disadvantages, particularly in terms of their conservation properties. For central fluxes we use averages of values just across the interface, i.e., and , where represents the averaging operator,
[TABLE]
for any function . On the other hand, using upwind fluxes requires solving a Riemann problem in a coordinate system local to that face. Consider a local coordinate system on the configuration space cell face, i.e., on . Here, is a unit vector normal to , and and are tangent vectors such that . Let and be electric and magnetic fields in this coordinate system. Then, assuming variations only along direction , Maxwell’s equations reduce to , , and the following uncoupled set of two equations for the tangential field components
[TABLE]
and
[TABLE]
Multiplying the first of each pair by and adding and subtracting from the second of that pair we obtain a set of four uncoupled advection equations
[TABLE]
and
[TABLE]
Hence, the solution to the Riemann problem with initial conditions
[TABLE]
for , and
[TABLE]
for at is simply
[TABLE]
and
[TABLE]
Rearranging these expressions shows that the upwind fields in the local face coordinate system are
[TABLE]
and
[TABLE]
where is the jump operator
[TABLE]
for any function . The solutions to the Riemann problem given by Eqns. (51)-(54) are identical to those presented in previous studies of Maxwell’s equations [37].
Remark 1**.**
Before proceeding to the conservation properties of our semi-discrete system, we note that the discretization of Maxwell’s equations given by Eqns. (36) and (37) does not include the constraints given by Eqns. (11) and (12), i.e., the divergence constraints in Maxwell’s equations. Thus, our algorithm may violate these constraints over the course of the simulation. Where appropriate in section 5, we will discuss how the violation of the divergence constraints in Maxwell’s equations manifests.
3.3 Conservation properties of the semi-discrete system
We proceed as we did with the continuous system, first considering whether the discrete system conserves number density.
Proposition 6**.**
The discrete scheme conserves total number of particles.
Proof.
Choosing in the discrete weak-form, Eq. (30), and summing over all phase-space cells , gives the desired conservation law since the volume integral vanishes for and the surface integrals are symmetric about cell interfaces. ∎
Before we move on to the norm, we consider the following Lemma,
Lemma 1**.**
Phase space incompressibility holds for the discrete system, i.e.,
[TABLE]
Proof.
For the specific discrete phase space flow in the Vlasov-Maxwell system, . Within a cell, Eq. (56) is zero since, as with the continuous system, has no configuration space dependence, and has no divergence in velocity space. The question is whether the jumps in across cell interfaces in phase space are accounted for by the scheme. Integrating Eq. (56) over a phase space cell , employing Gauss’ Law, and summing over cells,
[TABLE]
This result follows for the simple reason that the phase space flow is in fact continuous with respect to the surfaces considered, allowing us to pairwise cancel the integrand upon summation. For example, consider the configuration space component of the flow , . The velocity, , is continuous across configuration space surfaces because has no configuration space dependence. Likewise, the velocity space component of , is continuous across velocity space surfaces because and have no velocity space dependence, and in the term is the velocity coordinate, and thus is continuous. We note that this proof is specific to the phase space flow for the Vlasov-Maxwell system and in general may not hold for all systems. ∎
Proposition 7**.**
The discrete scheme decays the norm of the distribution function monotonically.
Proof.
Since the distribution function itself lies in the test space, we can set in Eq. (30). We then have,
[TABLE]
where we have used Lemma 1 to rewrite since phase space is incompressible, even in our discrete system. Now, summing over all cells and inserting our numerical flux, with a bit of algebra, we obtain,
[TABLE]
where the is determined by the value of in the Lax-Friedrichs flux. Note that the sum over surfaces is a sum over pair-wise surfaces, and we have included the contributions on both the left and right side of the surface. Recall that for the proof of phase space incompressibility for the discrete system, we exploited the fact that the phase space flow is continuous with respect to the surfaces considered, i.e., is continuous across configuration space surfaces and is continuous across velocity space surfaces. We can use this fact again to simplify our expression,
[TABLE]
So, the time evolution of the norm of the distribution function is bounded from above by a negative-definite expression. Thus, the norm of the distribution function is a monotonically decreasing function, and our scheme is stable. ∎
Corollary 1**.**
If the discrete distribution function remains positive definite, then the discrete scheme grows the discrete entropy monotonically444The monotonic growth of the discrete entropy is due to our convention in the definition of the entropy. If one drops the minus sign in the definition of the entropy, then the discrete entropy is a monotonically decreasing function, instead of an increasing function, if the discrete distribution function remains positive definite.,
[TABLE]
Proof.
Using the well known bound,
[TABLE]
we can see that , so long as remains a positive definite quantity, and thus is well-defined. Multiplying by then gives us the inequality,
[TABLE]
But the left-hand side is just the discrete entropy. Integrating over a phase space cell , summing over cells, and taking the time-derivative of both sides gives us an expression for the time evolution of the discrete entropy in our scheme,
[TABLE]
Now, we note that in Proposition 7 we have already proved that the norm of the discrete distribution function is a monotonically decaying function, and thus the negative of the norm is a monotonically increasing function, and by Proposition 6, the discrete system conserves particles. Thus, the discrete entropy is a monotonically increasing function. The significance of this is simply that, while we will demonstrate shortly that the semi-discrete scheme derived here conserves energy exactly, the scheme also has a small amount of “regularization” associated with it analogous to numerical diffusion present in the discretization of other equation systems via similar methods. Although, in this case, the regularization occurs in the complete phase space as opposed to the numerical diffusion present in discretizations of fluid equations. ∎
We will show below that the discrete DG scheme also conserves total energy exactly if central fluxes are used for Maxwell’s equations. The energy conservation is independent of the flux used in the Vlasov equations. We will first consider the conservation of energy for our discretization of Maxwell’s equations.
Lemma 2**.**
The semi-discrete scheme for Maxwell’s equations conserves electromagnetic energy exactly when using central fluxes, and is bounded when using upwind fluxes,
[TABLE]
By bounded, we mean that when using upwind fluxes, when the right hand side is positive, the electromagnetic energy will increase less than , and when the right hand side is negative the electromagnetic energy will decay more than .
Proof.
From the discrete weak-form of Maxwell’s equations, we need to compute equations for and . Since each component of the field lies in the selected test space, we take the -component of Eq. (36) and use as a test function, e.g., choose . Summing these three equations will give us an expression for the time-derivative of . We follow the same procedure for Eq. (37), which gives an expression for the time-derivative of . With a bit of algebra, we obtain,
[TABLE]
and
[TABLE]
We now multiply both equations by and add them. Since , we can combine the third terms of Eqns. (66) and (67)
[TABLE]
In the above result, note that upon integration by parts, we must use the field just inside the face. Hence, the evolution of the electromagnetic energy in a single cell becomes
[TABLE]
Exact Energy Conservation With Central Flux. Using central-fluxes to determine the interface fields, i.e., setting and , gives us,
[TABLE]
Upon summing over all configuration space cells, we see that the surface term vanishes as it has opposite signs for the two cells sharing an interface, leading to the desired discrete electromagnetic energy conservation equation,
[TABLE]
Energy Bound With Upwind Flux. To see what happens when using upwind fluxes, we transform the fields appearing in surface integral into the coordinate system. We can then write the third term in Eq. (69) as Using Eqns. (51)-(54) for the interface fields and summing over all configuration space cells, we then obtain
[TABLE]
Note that due to the symmetry of the terms, the central flux terms in Eqns. (51)-(54) have vanished on summing over all cells. Now consider the contribution of the term to the two cells adjoining some face. This term will be and . On summing over the two cells, this contribution will become . Likewise for the other electric and magnetic field coordinates. Hence, the surface terms, on summation, contribute non-positive quantities to the right-hand side, implying that
[TABLE]
∎
Lemma 3**.**
If belongs to the approximation space , then the semi-discrete scheme satisfies
[TABLE]
Note that the species index is implied, the sum over in the first term is over all phase space cells, and the sum over in the second term is over all configuration space cells.
Proof.
As we have assumed that , we can set in Eq. (30). This assumption gives us,
[TABLE]
where we have used the continuity of in the second term, so there is no need to evaluate the basis function just inside the cell. Summing over all species and all phase-space cells, the surface term vanishes as the interface flux contributes terms with opposite signs exactly like our density conservation proof, and we obtain the required identity
[TABLE]
Note that we have carried out the sum in velocity space to compute the current density, leaving an integration and sum over configuration space cells. This identity would hold if we also reduced the Lax-Friedrichs flux to either a central flux or upwind flux numerical flux function. We emphasize that this proof only requires because that is all that is required to substitute for the test function . This may appear deceptive as the continuous proof in Proposition 5 appears to contain terms which are higher order than . However, we note that the higher order terms in the continuous proof come from the substitution of the explicit expressions for , the phase space flow, whereas in the discrete proof presented here, we have left the discrete phase space flow as is to stress the fact that has its own basis function expansion. Thus, the higher order terms which are explicit in the continuous energy conservation proof are implicit here in the discrete energy conservation proof. If , then as well, and terms in the discrete phase space flow such as , the configuration space component of the phase space flow, can be exactly represented in terms of our basis function expansion. ∎
With these lemmas, the proof of total energy conservation is straightforward.
Proposition 8**.**
If central-fluxes are used for Maxwell equations, and if , the semi-discrete scheme conserves total (particles plus field) energy exactly.
Proof.
Using Lemma 2 and 3 demonstrates that the semi-discrete scheme conserves total energy exactly,
[TABLE]
∎
We remark that if upwind fluxes are used for Maxwell’s equations, then the energy will decay monotonically. This is because the bounding of the electromagnetic energy when using upwind fluxes results in a small amount of energy being lost from the system, i.e., when is positive, more energy is lost from the electromagnetic fields than goes into the particles, and when is negative, less energy goes into the electromagnetic fields than was given by the particles. Either way, the conversion of electromagnetic energy from, and to, particle energy results in some decrease of the total energy when upwinding is employed with Maxwell’s equations. We note though that other authors have demonstrated this loss of energy is small for higher order schemes such as the DG method employed here [38], and we confirm this fact by studying the energetics of many of the benchmarks presented in Section 5. Further, just as with Lemma 3, Proposition 8 requires that . By employing at least piecewise quadratic polynomials in velocity space, we can insure that is in our space of test functions.
Remark 2**.**
We note that while employing at least piecewise quadratic polynomials is a way to guarantee , we can also project onto piecewise linear basis functions. In that case, it will be a projected energy which is a conserved quantity and thus the semi-discrete system can conserve energy even when utilizing piecewise linear basis functions.
Lack Of Momentum Conservation. Note that we have made no mention of whether the discrete system conserves momentum. This omission is because the discrete system does not in fact conserve momentum. We can show momentum non-conservation by choosing and proceeding as we did with the continuous system,
[TABLE]
Since is continuous, upon summation over all phase space cells and species, we obtain,
[TABLE]
We can proceed exactly as we did with the continuous proof, but we note a key subtlety,
[TABLE]
Since the electric and magnetic fields are discontinuous across configuration space cell interfaces, we cannot use integration by parts to eliminate the latter two terms. In other words, integration by parts holds only locally and not over the whole domain due to the jumps in the fields across surfaces. However, it is important to note that from the form of this equation that in fact momentum conservation depends only weakly on velocity space resolution, and since the size of the discontinuities in the electric and magnetic fields decrease with increasing configuration space resolution, we can more strongly conserve momentum by increasing configuration space resolution.
3.4 The time discretization of the semi-discrete Vlasov-Maxwell system
Now, given the semi-discrete Vlasov-Maxwell system, all that remains is to specify a means of solving the resulting system of ordinary differential equations to advance the solution in time. We have implemented a third order strong-stability Runge-Kutta (SSP-RK3) method to generate a fully explicit update [39]. We find for a multi-stage method, SSP-RK3 is a good optimization of memory footprint considerations and accuracy for our time-discretization. For higher order time discretizations, the likely avenue is not more stages in a Runge-Kutta type scheme, but high order single step methods such as an Arbitrary DERivative (ADER) scheme like those employed in Balsara et al. [40].
Since the update is fully explicit, it is worth briefly mentioning what constraints on our time step exist, i.e., what Courant-Friedrich-Lewy (CFL) conditions exist for this system of equations. Given that we are solving Maxwell’s equations, a constraint in configuration space is the speed of light. In every problem of interest the speed of light will be the largest velocity in the system. This constraint can be written as,
[TABLE]
where is the speed of light, is the number of configuration space dimensions, and is the polynomial order of our expansion. This particular CFL condition is similar to the constraint for the finite-difference-time-domain (FDTD) discretization of Maxwell’s equations [41], but with an additional safety factor for stability which depends upon the number of configuration space dimensions and polynomial order of our basis expansion. We note that because we will always choose the speed of light to be the largest velocity in the system, Eq. (81) will be a more severe time step restriction than the constraint from the configuration space component of the phase space flow, i.e., the term. However, Eq. (81) is not the only CFL condition in a fully explicit update of the Vlasov-Maxwell system. Since we are solving for the distribution function on a combined configuration and velocity space grid, there is a corresponding maximum acceleration for a parcel of distribution function in velocity space. The CFL condition is then given by the phase space flow in velocity space,
[TABLE]
where the max on the phase space flow is to determine in which direction in velocity space is the acceleration the largest, and this time is the number of velocity space dimensions. In other words, the maximum component of the total Lorentz force, , will also set a constraint on the time step. It is not uncommon, depending upon simulation parameters, for this constraint to be more restrictive than the speed of light CFL. For example, the force can be quite large at the edge of velocity space for moderate magnetic field amplitudes depending upon how large the velocity grid is.
4 Efficient Implementation of Kinetic Update
Since the equation system we are solving is intrinsically high-dimensional, with many problems of interest requiring the solution of a four, five, or six dimensional hyperbolic equation, plus time, special care must be taken when dealing with the so-called “curse of dimensionality,” which refers to the exponential scaling of solver cost with dimension size. In a finite element framework, the “curse of dimensionality” is partly due to the polynomial expansion within a cell. A typical polynomial space employed is the Lagrange Tensor product space, formed by taking tensor products of a one-dimensional monomial expansion. Due to the nature of the tensor product, the number of degrees of freedom within a cell will scale as , where is the polynomial order and is the dimension. For four, five, and six dimensional simulations, this scaling quickly becomes intractable. As mentioned previously, we employ a reduced basis set commonly referred to in the literature as the Serendipity basis set, which retains the formal convergence order of the Lagrange Tensor basis, while eliminating many monomials from the tensor product expansion in higher dimensions.555It should be noted that, while this is true for structured grids, the convergence order of the Serendipity expansion on an unstructured grid is more subtle. This is to say, arbitrary refinements of an unstructured grid will destroy the convergence order of the Serendipity expansion, but with care the convergence order of the space can be maintained [42]. All work considered herein is on structured quadrilaterals. As such, the Serendipity basis’ scaling in higher dimensions is much less severe than the Lagrange Tensor basis, scaling like
[TABLE]
where is the number of polynomials in the expansion within a cell, and as before, is the polynomial order and is the dimension. The essential feature of the Serendipity polynomial space is that the monomials with “super-linear” degree greater than are not needed to maintain the formal convergence order of the tensor product expansion. By super-linear degree, we mean the degree of the monomial without counting the linear monomials. For example, in three dimensions, with polynomial order two, the Lagrange Tensor basis will have monomials of the form , but this monomial has superlinear degree four and is subsequently dropped by the Serendipity basis expansion. A comparison of the number of degrees of freedom within a cell for the two basis sets, given in Tables 1 and 2, reveals a stark contrast, especially for five and six dimensions even for relatively low polynomial order.
The actual construction of the algorithm now proceeds in two steps. Having chosen our polynomial space, we must choose a criteria for the creation of the polynomial expansion. Two common choices are a modal and nodal expansion. These two formulations can be shown to be equivalent mathematically, i.e., there exists a linear transformation between the two choices of basis expansions. However, their comparative computational construction is different. A modal expansion specifies only a set of polynomials, which can be as simple as for polynomial order , to more complicated special basis sets such as Legendre or Chebyshev polynomials. Here we choose the nodal expansion, specifying a set of nodes on a reference quadrilateral element for which a polynomial takes the value of 1 at one node, and 0 at all other nodes. The node location is somewhat arbitrary, though certain nodal configurations provide key computational advantages. Here we choose such a nodal layout, schematically given in 1D, 2D, and 3D in Figure 1.
This particular nodal layout is constructed such that every higher dimensional reference quadrilateral element is built from the lower dimensional reference quadrilateral elements, so that the lower dimensional faces of a reference quadrilateral element also form a unisolvent expansion, i.e., the polynomials local to the face form a complete basis of the solution space. For example, consider the pictorial representation of the 3D reference element. Each 2D face of the reference 3D element is exactly the 2D reference element for that particular polynomial order. This same recursive approach can be applied to higher dimensions as well666This fact is true in general, but higher polynomial orders may modify the lower dimensional reference elements such that the recursive algorithm is not quite as obvious as the one presented here. Just as polynomial order 4 introduces an interior node to a reference 2D element, so can higher polynomial orders introduce interior nodes to higher dimensional reference elements which would have to be taken into account in the recursive generation of the reference element. For up to polynomial order 4 though, every higher dimensional object can easily be generated as described, with 2D reference elements making up the faces of a 3D reference element, 3D reference elements making up the faces of a 4D reference element, and so on. We emphasize that this is likely sufficient for our problems since the Serendipity basis in four, five, and six dimensions, with polynomial order 4, involves the solution of a large number of degrees of freedom per cell, and due to the same performance and cost considerations that motivated the use of the Serendipity basis, we seek to avoid evolving thousands of degrees of freedom per cell., with the reference 4D element being comprised of reference 3D elements for each of the 4D element’s eight 3D faces, and a reference 5D element consisting of a reference 4D element for all ten 4D faces of a 5D element. This approach has the advantage of greatly simplifying surface integral calculations. Since every higher dimensional element is recursively generated from lower dimensional elements, every face of a higher dimensional element, the 2D faces in 3D, or the 4D faces in 5D, forms a unisolvent expansion for that surface. We thus only require the nodal information local to that face and can reduce the number of multiplications in the evaluation of the surface integrals by a somewhat sizable fraction. In 5D for instance, to advance the solution of the distribution function in time, one performs one 5D volume integral and ten 4D surface integrals. By requiring only the nodal information on a 4D surface, the number of degrees of freedom in evaluating the surface integral given in Table 2 is reduced by a factor of two to three.
The second step is then, given the choice of polynomial expansion, determining a means of evaluating the volume and surface integrals in the DG spatial discretization. Gaussian quadrature is a natural choice. The standard approach is to establish a tensor product of quadrature points in each direction for every dimension one wishes to integrate, and this approach will integrate exactly monomials of a particular order, for Gauss-Legendre or for Gauss-Lobatto for example, regardless of the dimension in which monomial varies. Here, we employ a Gauss-Legendre integration scheme, but we note that this strategy quickly becomes untenable for the same reason the Lagrange tensor polynomial space is prohibitively expensive for solving the Vlasov equation: the “curse of dimensionality.” For example, consider integrating the volume term in five dimensions with second order polynomials. Naively, one expects this to require the integration of monomials with degree in each dimension, because both and have polynomial expansions, thus requiring 4 quadrature points in each dimension, or a total of quadrature points, to avoid under-integrating the volume term in a cell. Given that the scaling of the computation of the volume integral in a cell is , the number of operations per phase space cell becomes quite large for modest polynomial orders in high dimensions. One way to reduce the number of operations dramatically would be to accept aliasing errors from under-integrating the volume and surface integrals in the DG discretization of the Vlasov equation. Accepting aliasing errors is a common approach to reducing the cost of fluid codes, where the nodes which specify the polynomial expansion in a cell will be co-located with quadrature points, and the interpolation to quadrature points operation is avoided, usually at the cost of under-integrating nonlinear terms. However, this is not an acceptable solution for numerically solving the Vlasov equation because the conservation relations derived in Section 3 rely on the integrals being computed correctly, i.e., the conservation relations in the Vlasov-Maxwell system are implicit, as opposed to the explicit conservation relations which form fluid systems of equations. Under-integration of the DG discretization of the Vlasov equation can lead to poor particle and energy conservation, which can eventually lead to numerical instability.
Fortunately, the naive requirement to avoid under-integration of terms in our DG discretization is quite far from reality due to the structure present in the Vlasov equation. Consider the advection in velocity space,
[TABLE]
To avoid over-complicating the notation, we have dropped the species index. We can see clearly that in fact, integrating monomials of degree is only required in configuration space, and in velocity space, we require at most integrating monomials with degree . We can thus employ an anisotropic grid for performing the quadrature of the advection in velocity space, using only the minimum number of quadrature points required along each direction of integration. Table 3 considers the impact anisotropic quadrature has on a few combinations of velocity space and configuration space dimensions. While there is no gain for polynomial order 1, the gain is a moderate improvement for other combinations compared to isotropic quadrature, as shown in Table 4. A similar reduction in the number of quadrature points required can be demonstrated for the surface integrals. There is an even larger reduction in the number of operations with the configuration space component of the advection in phase space on a structured, Cartesian, grid. Since is a coordinate, it requires no expansion in basis functions.
Upon rewriting the term then as,
[TABLE]
we note that both of these terms are such that the integrals can be pre-computed. Here, is velocity in the center of the cell. In other words, upon expansion of the distribution function in basis functions given in Eq. (31), to update the distribution function in configuration space, we only require two matrices,
[TABLE]
along with the cell center coordinate. We remark that is commonly referred to as the grad-stiffness matrix in finite element literature, weighted by the constant cell center velocity coordinate. In rewriting the configuration space update in this way we reduce the update from operations to operations, where is the number of quadrature points required to exactly integrate either the volume or surface configuration space terms. for the configuration space update scales like , so the gain in writing the update this way is identical to the gain we obtained from employing the Serendipity polynomial space instead of the Lagrange Tensor basis. We note that similar manipulations can be done for the surface integrals in the DG discretization of the Vlasov equation to reduce the number of operations in those components of the update. For example, the quadrature required for the force term surface integrals is reduced not only by the exploitation of anisotropic quadrature requirements, but also by the fact that since as mentioned previously, each set of surface nodes is itself a unisolvent expansion on that surface, we use only that nodal information in the evaluation of the surface integrals. Likewise, the streaming term surface integrals can be precomputed in a similar fashion to the volume integrals due to the fact that the configuration space component of the Lax-Friedrichs flux can be computed analytically,
[TABLE]
For the Cartesian grids employed here is simple to calculate, e.g., surface integrals computed at constant in configuration space just involve the velocity coordinate . Given the analytic form of the Lax-Friedrichs flux, similar matrices to Eq. (86) and Eq. (87) can be derived and multiplied by or depending upon the sign of the velocity,
[TABLE]
where are the surface basis functions, i.e., basis functions local to a surface defined by the nodal information on a surface, and we have introduced the surface element to emphasize these surface integrals are for evaluating the configuration space component of the phase space flow and thus surfaces at some constant configuration space dimension. Note that is the basis function of the full volume expansion, but evaluated at the surface of interest so that the dimensions of this matrix are . While individually, each of these reductions in the number of operations per time step is modest, taken altogether, they provide a substantial gain over a traditional finite element approach using a tensor product basis.
A final note on the implementation involves the coupling of the Vlasov solver and the solver for Maxwell’s equations. Since we require the first velocity moment to advance Maxwell’s equations in time, it is worthwhile to consider efficient algorithms for the calculation of velocity moments. For simplicity, let us first consider how one would calculate the zeroth velocity moment, the density. The density in this case is related to the charge density defined earlier,
[TABLE]
In the discrete system, the density belongs to the same space as the electromagnetic fields , and thus the expansion of the density involves the basis expansion in only configuration space, i.e.,
[TABLE]
The goal is then to determine . We do this by multiplying the definition of the discrete zeroth moment on both sides by and integrating over configuration space. Note that by discrete zeroth moment, we mean that since phase space has been partitioned into a phase space mesh with cells , Eq. (91), which is a global operation in velocity space, much also be partitioned. We do this by summing over the contributions of the velocity integral inside each cell , but only for a fixed position in configuration space, i.e.,
[TABLE]
where our choice of notation implies we will sum only over the space orthogonal to configuration space for each cell . On the left hand side of the Eq. (94), we note that the integral is what is commonly referred to in finite element literature as the mass matrix, but only in configuration space,
[TABLE]
The right hand side is what we define as the zeroth moment matrix,
[TABLE]
We remark that, on uniform, structured grids, the zeroth moment matrix is the same in every cell and thus the algorithm defined above can be pre-computed with minimal memory storage. Non-uniform structured grids would require a small generalization in the form of a constant scale factor. Higher moments, such as the first velocity moment needed to compute the current for Maxwell’s equations,
[TABLE]
are slightly more subtle. Note that a first moment matrix defined as,
[TABLE]
would be different in every velocity space cell, even on uniform, structured grids. A similar rearrangement to the one employed for the configuration space component of the advection in phase space provides a significant savings in terms of memory,
[TABLE]
On a Cartesian mesh, it is most obvious how this manifests algorithmically if we break up the update into its vector components. In other words, if the system has three velocity dimensions, then
[TABLE]
where . This result can be generalized to higher moments such as the stress tensor and energy density,
[TABLE]
and beyond.
Before moving on to benchmarks, we summarize the results of this implementation section by synthesizing all the steps above into a coherent algorithm, i.e., how each operation described manifests in a single Runge-Kutta stage. Define and to be the values of the distribution function before and after the Runge-Kutta stage respectively.
Compute current from old value of the distribution function
- (a)
The current is computed by summing over the contribution from the first moment for each species,
[TABLE]
where is the inverse of the matrix given by Eq. (95) and is given by Eq. (100). 2. 2.
Advance the Vlasov equation one time-step for each species.
- (a)
Compute the volume integral component for the streaming term,
[TABLE]
where is given by Eq. (86) and is given by Eq. (87). The matrix in this case is related to Eq. (95), but this matrix involves phase space basis functions,
[TABLE]
- (b)
To obtain the volume integral component of the force term, interpolate the old distribution function, , and the old velocity component of the phase space flow, , onto Gauss-Legendre abscissas, but with the specified anisotropy given by Table 3, and evaluate the integrals using Gauss-Legendre quadrature.
- (c)
Increment with this additional contribution multiplied by .
- (d)
Compute the surface integral component for the streaming term,
[TABLE]
where is given by Eq. (89) and is given by Eq. (90). The matrix is identical to the matrix used in the volume component of the streaming term, i.e., the inverse of Eq. (105).
- (e)
To obtain the surface integral component of the force term, we consider the nodal components of and the velocity component of the phase space flow, , local to the surface of interest and interpolate that nodal information onto Gauss-Legendre abscissas, again with the specified anisotropy so that the minimum number of quadrature points required are used along configuration and velocity space dimensions. We then evaluate the Lax-Friedrichs flux at each quadrature point, taking into account the sign of so that the upwind direction is computed correctly, and compute the surface integrals using Gauss-Legendre quadrature.
- (f)
Finally, as before, we increment with this result multiplied by . 3. 3.
Advance Maxwell’s Equations one time-step.
- (a)
Since Maxwell’s equations are linear, we employ standard techniques for linear advection equations, with the caveat that the form of the numerical flux is either central fluxes, e.g. the simple average of the left and right fields, or upwind, as defined by Eqns. (51-54).
- (b)
We then increment the component of the electric field with the current derived from given by Eq. (103), i.e.,
[TABLE]
This procedure is repeated for each Runge-Kutta stage and is extensible to an arbitrary order multi-step method, though the memory usage of a multi-stage method beyond SSP-RK3 may be prohibitive. We note that the inverses of the various mass matrices which appear do not need to be multiplied through for each time-step on the structured, Cartesian, grids considered in this paper and can be pre-factored into the form of the update step. Further, it is worth acknowledging why for the subset of meshes considered herein we do not pre-evaluate all the integrals present, since polynomials on structured grids can be pre-computed easily. The simple answer is two-fold. The first component of the answer is that for the nodal basis expansions considered, the pre-computation of all the integrals does not provide significant savings in terms of the number of operations required to perform a Runge-Kutta stage. Roughly, this algorithms most expensive component, the non-linear force term, scales like where is the number of quadrature points, and is the number of basis functions in the volume expansion. Pre-computation of this component of the update reduces the number of multiplications required to where is the number of basis functions in configuration space, i.e., the number of basis functions employed in the expansion of Maxwell’s equations and the moments of the distribution function. Upon comparison of these two scalings, we find at best for the high dimensional problems considered a factor of 2 savings in the number of operations. However, the second component of the answer is that the memory requirements of the pre-computation update are more severe due to the need to store a large rank three tensor for each velocity dimension. For the form of the basis expansion we employ in this paper, we found that the additional memory requirements were not worth the small reduction in the number of multiplications. This may not be true in general, as a modal basis expansion in orthogonal polynomials could have a significant reduction in the memory footprint due to the resulting sparsity of the tensor. We will consider extensions of the algorithm to different basis expansions in a subsequent publication.
5 Applications
In the following sections we present a number of benchmarks and studies as verification of our discretization of the Vlasov-Maxwell system. The following set of simulations is by no means exhaustive. Further verification of the Vlasov-Maxwell solver can be found in [23, 24], which employ the Vlasov-Maxwell solver to study the plasma sheath and the current filamentation instability (CFI), a drift induced instability caused by counter-streaming beams, which can serve to amplify magnetic fields. We note in particular the agreement [23] finds between linear theory of the CFI for long wavelengths and the linear growth of the instability in simulations, as well as the comparable nonlinear saturation amplitudes between two-fluid and Vlasov simulations of the instability. Further details of the nonlinear regime of the CFI, including a detailed analysis of the nonlinear saturation of the CFI can be found in [24].
5.1 Scaling studies
We briefly note a few scaling benchmarks of the Vlasov-Maxwell solver. Using a 1X3V setup, meaning one configuration space dimension and three velocity dimensions , with cells and polynomial order 4 with the Serendipity space, 136 degrees of freedom per cell, a strong scaling study demonstrates the solver’s ability to scale to thousands of processors. We can likewise perform a weak scaling study, starting from a box size of box on 16 processors and increasing to a the box size of on 2048 processors. Each simulation is run for hundreds of time steps to minimize the ratio of time spent initializing the simulation to time spent taking time steps. Results for these two scaling studies are given in Figure 2.
5.2 Electrostatic applications
Our first tests of the algorithm will be on systems for which the magnetic field is not dynamically important. While we solve Maxwell’s equations, Eqns. (9)-(12), in their entirety, it is not incorrect to draw analogy between the tests presented below, and the so-called Vlasov-Ampere system, where Maxwell’s equations reduce to,
[TABLE]
For all the following benchmarks, the distribution function is initialized as a drifting Maxwellian,
[TABLE]
where is the number density defined in Eq. (91), and are the flow and thermal velocity of species respectively. The flow naturally follows from Eq. (97),
[TABLE]
and the thermal velocity can be compute from the second moment of the distribution function,
[TABLE]
where the constant is determined by the number of velocity dimensions of the simulation, either one, two, or three. Alternatively, the thermal velocity can be defined in terms of a specified temperature,
[TABLE]
where is the temperature of the species.
5.2.1 Momentum and energy conservation tests
We first consider a simple problem to demonstrate numerically the conservation laws derived earlier. We will initialize a distribution function with strong asymmetric flows to determine how well momentum and energy are conserved in our implementation of the energy, but not momentum, conserving DG algorithm. In one configuration space dimension () and one velocity space dimension (), or 1X1V, we initialize a Maxwellian distribution function with density
[TABLE]
for both the protons and electrons. The length of the box is , where is the proton plasma frequency. We choose , , and , where is the electron inertial length. The drift is given by for both species. Other relevant parameters are , , and . The velocity space extents for the electrons are and the proton velocity space extents are . All simulations are run with . A number of CFL conditions and configuration space resolutions are employed to determine the convergence order of our time-stepping scheme and the convergence of the momentum errors. The boundary conditions in configuration space are periodic, while the boundary conditions in velocity space are zero flux to guarantee that there are no boundary losses for these tests of energy and momentum conservation. While it may appear that the density profile in Eq. (113) is not periodic, the values of and localize the density profile around since the profile varies on electron scales, while the domain size is with respect to proton scales. This scale separation reduces the violation of periodic boundary conditions to a vanishingly small number, . The simulations are run to where is the electron plasma frequency. This length of time corresponds to approximately ten thousand time steps for the runs with the largest time steps. Results are plotted in Figure 3. Since we are either halving the time-step or doubling the resolution, we define the order of convergence as
[TABLE]
where is the relative change of a quantity, either momentum or energy, in a particular simulation. In a set of simulations which are “converging,” so the convergence order is a positive number. We note that the energy conservation result can be obtained by either decreasing the size of the time step in an simulation, as shown, or by refining the grid. The energy errors for are roughly equivalent to the energy errors shown for the smaller time-step. For the energy conservation test, we find the orders of convergence are 2.12, 2.70, and 2.75 for each successive halving of the CFL number for an simulation. These orders of convergence are comparable with the expected third order scaling from our third order strong-stability-preserving Runge-Kutta method. For the momentum conservation test, we find the convergence orders of the polynomial order 2 simulations are 1.46, 2.27, and 3.76, and the convergence orders of the polynomial order 3 simulations are 1.99, 3.08, and 6.04, for a refining of the grid from to cells in configuration space. We can thus say the momentum conservation tests demonstrates super-linear scaling, i.e., greater than , for continual refinement of the grid.
We can also use this same initial condition to examine the behavior of quantities we expect to decay, or grow, throughout the simulation, such as the norm of the distribution function and the discrete entropy . However, we should point out that a crucial assumption in our proof of the growth of the discrete entropy was that the distribution function remain positive definite. Due to the nature of our DG scheme, this cannot be guaranteed. In fact, for highly resolved grids, while the positivity violations decrease in magnitude, they cannot be completely eliminated. The reason for this is the magnitude of the distribution function is so small at the edges of velocity space, that the generation of round-off error positivity violations is inevitable. Fortunately, as can be seen in Figure 3, small positivity violations do not affect the stability or robustness of our algorithm. This is not surprising, as even moderate positivity violations are tolerable so long as the moments of the distribution function remain positive. It would take substantial negativity in the distribution function to drive moments such as the density or energy negative and create the same numerical stability issues often found in fluid codes. Nonetheless, since we would like to at least comment on a quantity such as the discrete entropy even with the presence of positivity violations, we have run the same initial conditions with more resolution and defined an “entropy-like” quantity, , which we expect will behave similarly to the discrete entropy in the limit that positivity violations go to zero. We have plotted the norm of the distribution function and integrated over all of phase space for both protons and electrons in Figure 4 for highly resolved grids, with polynomial order 2, demonstrating, at least approximately, the behavior we expect based on Proposition 7 and Corollary 1.
Finally, we reiterate that our discretization of Maxwell’s equations is not guaranteed to preserve the divergence relations, in this particular case Eq. (11), and thus we would like to examine how well this simple initial condition preserves this constraint on the electric field. Since our system only has one dimension in configuration space, in a discrete sense we are concerned with how well the relation,
[TABLE]
is satisfied at each point in configuration space. We have plotted Eq. (115) in Figure 5, normalized to the maximum value of the divergence of the electric field, at for a number of resolutions in configuration space with polynomial order 2. We have also plotted the value of Eq. (115), integrated over all of configuration space, versus time, in Figure 5. The derivative of the electric field is computed using the analytic derivatives of the basis function expansion within a cell. The above normalization is somewhat arbitrary, but it does illustrate that the error in the divergence constraint is small compared to the largest variation in the electric field, and thus the errors in the divergence relation are sub-dominant to the dynamics. Further, it is a comfort to note that the errors decrease dramatically with increasing resolution, even given the complicated coupling at play since one of the sources of the divergence errors comes from the computation of the moments of the distribution function, and the charge density does not appear anywhere in the equation system we are explicitly evolving.
5.2.2 Landau damping of Langmuir waves
Consider a plasma, or Langmuir, wave propagating in a plasma of protons and electrons whose distribution functions are given by Maxwellians. Langmuir waves are dispersive waves, with a dispersion relation given by
[TABLE]
in the limit that the proton mass is much larger than the electron mass and the protons can thus be considered immobile. Here, is the electron thermal velocity, is the electron plasma frequency, and is the electron Debye length. is the plasma dispersion function, defined as
[TABLE]
with the derivative of the plasma dispersion function given by
[TABLE]
An application of complex integration techniques shows that depending on the sign of the largest imaginary component of the frequency , the wave is either unstable and will grow with time, or will damp away, a phenomenon known as Landau damping. For Langmuir waves propagating in a Maxwellian plasma of protons and electrons, the waves quickly damp away. Using a 1X1V setup, we can initialize Langmuir waves in the Vlasov-Maxwell system with a small density perturbation and the corresponding electric field to support this density perturbation,
[TABLE]
where , is the size of the perturbation, and is the wavenumber of the wave. The electric charge and permittivity of free space are included in the electric field to satisfy Eq. (11). Choosing allows us to compare with the analytic theory described above. The box size is set to so exactly one wavelength fits in the domain. Specific parameters for these runs are: , , , and . For the proton species, the velocity space extents are , and for the electrons, the velocity space extents are . The boundary conditions in configuration space are periodic, while the boundary conditions in velocity space are zero flux to preserve density and energy conservation in the semi-discrete system. The resolution is chosen for each simulation to adequately resolve the Debye length in configuration space and to mitigate numerical recurrence in velocity space. By numerical recurrence, we refer to the process by which the collisionless system artificially “un-mixes” if the distribution function forms structure at the velocity space grid scale. Numerical recurrence is inevitable with finite velocity resolution for this particular problem as the Landau damping of the wave will create smaller and smaller velocity space structure through the phase-mixing of the wave. We could completely eliminate this issue with a diffusive process in velocity space, such as a collision operator. Here, we choose ample velocity resolution so that the wave damps enough for us to extract a clean damping rate and frequency for the wave initialized. We find for the longest wavelengths, using polynomial order 2, a resolution of 64 points in configuration space adequately resolves the Debye length, and 128 points in velocity space permits the wave to phase-mix sufficiently to extract damping rates. The evolution of the electromagnetic energy, as well as the other components of the energy, in a prototypical simulation is given in Figure 6. Comparisons of a number of Vlasov-Maxwell simulations with theory for both the damping rates and the frequencies of the waves are given in Figure 7. For the theory, we solve Eq. (116) using a root-finding technique. We emphasize that we solve the Vlasov-Maxwell system in its entirety, including the nonlinear term, and for both a proton and electron species. With the above simulation parameters, the plasma waves damp entirely on the electron species, so the approximation that the protons are essentially immobile in our dispersion relation holds to a high precision. We also wish to note that the resolution of 64 points in configuration space is not required for every simulation. For example, the prototypical simulation presented in Figure 6 uses only 16 points in configuration space, or approximately one grid cell per Debye length. As long as the gradients are properly resolved, the Vlasov-Maxwell discretization is extremely robust. Subsequent simulations have a large amount of variation in their resolution compared to the Debye length. While Section 5.2.3 uses a similar resolution of approximately one grid cell per Debye length, a grid cell in the simulation presented in Section 5.3.2 contains 40 Debye lengths.
5.2.3 Electrostatic shocks
When two counter propagating, but otherwise identical, slabs of plasma collide at speeds greater than the proton sound speed , a potential well can be created, leading to trapping of electrons and the formation of shocks from nonlinear steepening of wave fronts. In the absence of collisions, shock formation only occurs up to a critical mach number , where is the drift velocity [43]. To demonstrate that the collisionless Vlasov-Maxwell system behaves as expected, we perform two simulations, one with and one with , and compare the results to analogous five moment two fluid simulations given by the algorithm described in Hakim 2006 [17]. The parameters of the two simulations are as follows: , , , and . The drifts of both the protons and electrons are equal, so as not to initialize any current. The velocity space extents of the electrons are , while the proton velocity space grid is , where is the magnitude of the specified drift. The proton velocity grid is made larger by the inclusion of the drift in its extents, but including the drifts ensures that the moments of the proton distribution function are calculated accurately. Including the drifts in the electron velocity space extents are unnecessary since the drifts are much smaller than the electron thermal velocity. The boundary conditions in velocity space are again zero flux, but in this case the configuration space boundary conditions are open to avoid any pollution of the region downstream of the shock due to waves launched by the initial condition propagating towards the wall. The resolution of the Vlasov-Maxwell simulations is with polynomial order 2. The left half of the domain is initialized with the distribution function propagating to the right, while the distribution function on the right propagates to the left. The drifting plasmas thus collide in the middle of the domain. At , we note the following comparisons between the two Vlasov runs with the two fluid simulations in Figure 8. Essential features of the shock for such as the mass density and momentum for both species compare well, while the particle energy, especially for the protons, shows poor agreement between the two models. In addition, while the fluid simulation still clearly shocks, the broadening of the densities in the Vlasov run demonstrates that no shock forms for this high of a Mach number in the collisionless system, in agreement with the aforementioned previous studies [43]. Plots of the distribution function in Figure 9 further illustrate some of the key differences between the fluid and kinetic models. For example, the protons in the are very far from thermal equilibrium, having experienced a significant amount of particle acceleration, thus easily explaining why the proton particle energy does not compare very well.
5.3 Electromagnetic applications
We now move on to more general applications, which involve both the electric and magnetic fields in Maxwell’s equations, Eqns. (9)-(12). For all benchmarks the distribution function is again initialized as a drifting Maxwellian as in Eq. (109).
5.3.1 Vlasov-Maxwell Riemann problem
A common robustness test of algorithms for fluid equations involves so-called “Riemann” problems, of which the Brio-Wu shock tube is a famous example [44]. Here we perform a 1X3V, , simulation of a Riemann problem with our Vlasov-Maxwell solver and compare with a case previously studied with a two-fluid algorithm. We have again employed the five moment two-fluid algorithm described in Hakim 2006 [17] to compare the fluid moments between the two simulations. The initial condition is given by
[TABLE]
where the left values are the fluid moments initialized on the left half of the domain, and the right values are the fluid moments initialized on the right half of the domain. We note that , the energy, is defined in Eq. (102). The full mass ratio, , is employed, and the domain size is where is the proton skin depth. All parameter definitions are with respect to values on the left half of the domain, i.e., is defined using the proton plasma frequency on the left half of the domain. The simulation is run for , where
[TABLE]
is the proton cyclotron frequency. We note that the cyclotron frequency is the same on both the left and right halves of the domain since the magnetic field has the same magnitude on each side of the interface. The velocity space extents of the Vlasov-Maxwell run are , where is defined by the higher temperature on the left side of the domain for each species, and zero-flux boundary conditions are employed in velocity space. As in the electrostatic shock problem, open boundary conditions are employed in configuration space. The speed of light is chosen such that , so that the edge of velocity space for the electrons is sub-luminal. The resolution of the Vlasov-Maxwell simulation is with polynomial order 2. We can compare some of the aforementioned fluid moments between the Vlasov-Maxwell and the same five moment two fluid model employed in Section 5.2.3. Results are shown in Figure 10.
We note some key differences between the two models due to the more complete description of the plasma the Vlasov-Maxwell system provides. The presence of a larger amount of structure in the two-fluid model is owed to the fact that in the fluid model, the dispersive waves launched by the sharp gradient in the density and the temperature can only propagate and nonlinearly interact with each other. However, these dispersive waves in the Vlasov-Maxwell system can also damp on the particles in the plasma, giving their energy to the plasma particles via a variety of wave-particle resonances, such as Landau damping. The presence of wave-particle resonance processes is likely the reason for the striking difference between the electron particle energy between the two models, where the profiles are similar in terms of global structure, but the Vlasov-Maxwell energy density profile is shifted up. We note that this shift is indeed due to the energization of the electrons and not due to any sort of energy conservation violations, as energy is conserved for this run to within 1 percent, with the only losses due to the fact that the configuration space boundary conditions are open. Further evidence for this damping can be seen by examining the electron distribution function. Isosurfaces of the electron distribution function are plotted in Figure 11. The ripples on the distribution function where the layer is initialized are a characteristic of significant phase-mixing, a key component of the damping of these waves. Evidence for phase-mixing is especially clear when considering 2D and 1D cuts of the distribution function, plotted in Figures 12 and 13.
5.3.2 The Orszag-Tang vortex
Our final benchmark is a 2X3V, , calculation of the Orszag-Tang vortex [45], a famous benchmark of magnetohydrodynamic (MHD) codes, and an extensively used simulation setup to study the physics of turbulence in 2D [46, 47]. The Orszag-Tang Vortex initial condition is given by:
[TABLE]
where we note the slight modification from the standard MHD initial condition is required in the equation for the electric field since we are now evolving a multi-species plasma. here is an out-of-plane guide field. The following parameters are employed: , where is the proton plasma beta, and is the Alfvén speed for species , . The amplitudes for the initial condition are, . is taken to be 1.0. The size of our box in configuration space is and the velocity space extents are . Periodic boundary conditions are employed in configuration space and zero flux boundary conditions are employed in velocity space. The resolution of the simulation is with polynomial order 2. Our resolution is such that, with polynomial order 2, we resolve the electron skin depth, . Schematically, the initial condition is given in Figure 14. For comparison, an analogous 10 moment two fluid calculation is run, though with larger resolution to avoid grid scale diffusion in the fluid algorithm. By 10 moments, we mean that the algorithm solves for all six independent components of the pressure tensor for both species. The algorithm description is given in Hakim 2008 [18]. After some time, the flow and magnetic field become turbulent due to the twisting of the field and flow. In Figure 15, the out-of-plane current and magnetic field contours are compared between the two models, and qualitative agreement is found.
Though many large scale features between the two models agree as they should, there are key physical differences between the two models which highlight the more complete description of the plasma the Vlasov-Maxwell system provides. Since the 10 moment two fluid model can be thought of as a truncation of the Vlasov-Maxwell model, it is apparent from the larger currents in the Vlasov-Maxwell simulation that there is physics missing from the fluid model, either due to the higher fluid moments beyond the pressure tensor, or the closure of the fluid system, which causes the thinning current sheets in certain regions of the domain. It is also instructive to consider the energetics and spectra of the two simulations, plotted in Figures 16 and 17 respectively.
Of course, one of the most interesting diagnostics remains: plots of the distribution function. Using a k-filter, where the distribution function is Fourier transformed, all Fourier coefficients except a single and are set to zero, and the distribution function is inverse Fourier transformed back to real space, one can consider how the distribution function is evolving at various scales. It is particularly noteworthy in Figure 18 that a large amount of structure is contained in the Vlasov-Maxwell system. Such features could easily be smoothed away in a standard PIC simulation due to noise and coarse graining.
A natural follow up question is then, what does such signal provide? Beyond the ability to more easily compare with observations and experiments if performing verification studies of the solver, having noise-free distribution functions increases the number of tools available for analysis of a simulation. For example, consider a diagnostic called a field-particle correlation [48, 49, 50]. There are three field-particle-correlations which we will consider here, defined as follows:
[TABLE]
where are directions with respect to the local magnetic field, is the correlation time, and is the configuration space point at which the correlation is being performed. The function is defined as,
[TABLE]
The correlation time determines the summation extents, i.e., , , while the index determines the specific time at which the correlation is being measured. In other words, we determine the correlation at time over a specified correlation time by partitioning the correlation time into chunks. The power of this diagnostic lies in its ability to ascertain secular energy exchange between the particles and fields, which may be occurring amidst complicated dynamics such as the broadband turbulence found in the Orszag-Tang set up. For example, the first correlation, Eq. (128), provides a proxy for determining whether Landau damping or parallel acceleration may be taking place at the specified point. Landau damping of a particular wave mode can be determined with a judicious choice of correlation time. The other two correlations, Eqns. (129)-(130), provide proxies for two other types of wave particle resonances: transit time damping and cyclotron damping. There are two ways to represent this data, as stack plots in , for analyzing how much parallel phase mixing may be occurring, or in at a specified time as a way of ascertaining where the wave-particle resonances may be occurring, and thus narrowing down what wave modes may be present in the simulation. These two ways are plotted in Figures 19 and 20. The level of signal is striking, as determining this amplitude of a correlation could be challenging in a standard particle-in-cell simulation without a sizable number of particles per cell, especially given the gradients in velocity space of the distribution function present in the field-particle-correlation definitions, which tend to exacerbate counting noise issues. A future publication will explore the significance of the correlations presented here.
We wish to make a couple final notes here about the Orszag-Tang simulation. One is that, out of all the benchmarks presented, it is perhaps most readily apparent that the small positivity violations described as inevitable in Section 5.2.1 are present in this simulation, e.g. see the upper panels of Figure 18. We reiterate that this behavior is, of course, unphysical. We also want to emphasize again that the fluid moments such as the density and energy remain positive definite for the entire simulation in spite of these positivity violations in the distribution function, thus further adding credibility to our claim that small positivity violations do not pose as severe of a numerical threat as positivity violations do in standard fluid codes. Nonetheless, the unphysical nature of a negative distribution function, if sufficiently large, is enough to give one pause about any analysis one wishes to perform. Unfortunately, as we move to higher and higher dimensions, increasing the velocity resolution dramatically to reduce these positivity violations becomes a more intractable solution. Limiters which preserve the positivity of particular quantities in fluid algorithms are incredibly common, but due to the fact that conservation relations are implicit in the Vlasov-Maxwell system, fixing positivity in a kinetic system has additional subtleties. In Cagas et al. [24], they employ the Moe-Rossmanith-Seal limiter [51] with the Vlasov-Maxwell solver,
[TABLE]
where is the cell average of the function , and is defined as,
[TABLE]
Such a limiter shifts the distribution function upward in a cell to guarantee positivity while keeping the cell average, and thus the number of particles, constant. However, energy conservation will be broken by this operation, and it could additionally make momentum non-conservation worse. In fact, Cagas et al. [24] notes that the operation of fixing positivity gives rise to one to two percent energy conservation errors due to the increase in the particle energy. Conserving momentum and energy after the application of such a limiter can be done by scaling the distribution function by an additional function
[TABLE]
where
[TABLE]
and is the second moment of the distribution function prior to the application of the limiter, while is the second moment of the rescaled distribution function after fixing positivity violations. Unfortunately, such an operation is nonlocal in velocity space, and it is as yet unclear what effect such an operation has on the numerical stability of the algorithm. A more natural limiter may involve modification of the differential operator, as opposed to the solution, to prevent the differential operator from ever generating unphysical positivity violations. By limiting the solution via changes to the differential operator, we could retain our conservation relations, but additional testing would be needed to determine if a modified operator still has the same convergence properties and produces the correct physical behavior of the Vlasov-Maxwell system. We defer a systematic study of limiters for kinetic systems to a future publication and for now comment that the physics of the Orszag-Tang simulation does not appears to be significantly impacted by the presence of positivity violations in the distribution function.
The second note is that, on the coarse resolution explored here, it is even more important to discuss the divergence errors which will arise from how we have discretized Maxwell’s equations. We focus in particular on the errors in the divergence of the magnetic field, since for the plasma parameters employed, the dynamics of the Orszag-Tang turbulence are magnetically dominated. In Figure 21, we plot the divergence of the magnetic field normalized to the maximum of the curl of the magnetic field. The choice of normalization is due to how important the curl of the magnetic field is for the development of the current layers seen in the evolution of the out-of-plane current plotted in Figure 15.
Since the errors in Figure 21 appear to be small, at least when compared to a quantity of dynamical interest, we thus expect that these divergence errors, like the small positivity violations, have not significantly impacted the physics present in this turbulence simulation. However, it is worth acknowledging how we might eliminate, or at least further suppress, these sorts of errors since 2 percent errors are not completely trivial. One method of mitigating divergence errors is through so-called divergence cleaning schemes, such as Maxwell’s equations in perfectly hyperbolic form [52, 53, 54]. This system of equations adds the two divergence relations as dynamical equations, with their time evolution represented by correction potential which contain the accumulated errors in the divergence of the electric and magnetic fields,
[TABLE]
where are dimensionless parameters which determine the speed at which errors propagate in the system. While this scheme is perfectly reasonable in vacuum, where the charge density, , is zero, for a plasma, the correction potentials may modify the total energy, and thus the conserved quantity, or worse, they may modify the Lorentz force in a non-trivial way. Alternatively, reconstruction methods have been quite fruitful in magnetohydrodynamics [55, 56, 57], where one only worries about the divergence of the magnetic field, and they have been recently extended to finite volume discretizations of relativistic two-fluid equations, as long as one explicitly evolves an equation for the charge density [40]. Finally, extensions to curl-conforming methods such as Nedelec finite elements, popular in the continuous Galerkin community, have been very recently achieved for both the discontinuous Galerkin method and equation systems with explicit sources, such as the Vlasov-Maxwell system [58, 59, 60]. These latter methods are particularly exciting, as they present a possible method of completely avoiding divergence errors, as opposed to just mitigating them, by virtue of the DG discretization containing a discrete version of the continuous vector identity , where is some vector field.
6 Summary and Future Work
In this paper, we have described a new algorithm for the numerical solution of plasma kinetic equations. We have demonstrated that the semi-discrete discontinuous Galerkin (DG) discretization of the Vlasov-Maxwell system conserves particles exactly and conserves energy exactly when central fluxes are used for Maxwell’s equations, independent of the choice of numerical flux function for the Vlasov equation. We have discussed at length an efficient implementation of the kinetic algorithm, which significantly reduces the number of multiplications required in a nodal basis expansion. This economic formulation of the algorithm scales incredibly well on standard computing architecture and has been used to verify the solver against common benchmarks, both electrostatic and fully electromagnetic, from only one configuration and velocity dimension, to a five dimensional calculation of turbulence.
There remain many avenues of future work to further improve the Vlasov-Maxwell solver. As discussed in Section 3.4, the kinetic equation’s CFL condition can be thought of as two distinct restrictions on the time step, a maximum velocity and a maximum acceleration, where the latter CFL can quite easily be more restrictive than the former. For instance, in the five dimensional turbulence run, because the magnetic field’s amplitude is so large, the acceleration CFL is as much as two times more restrictive than the speed of light CFL with the given parameters. This constraint on the time-step arises despite being significantly non-relativistic, since . Splitting the update into a configuration space and velocity space update, and treating just the velocity space update implicitly, would pay huge dividends. Further, as we hinted at in the summary of the algorithm, a reformulation of the DG update in terms of a modal update is an interesting path to pursue. While the nodal expansion often provides the simplest means of generalizing to complicated geometries and very sophisticated boundary conditions, for many astrophysical systems, simplicity is usually easily justified. As such, a modal expansion on structured, Cartesian grids, where a Gram-Schimdt orthogonalization procedure could be employed to produce an orthonormal basis expansion would produce a dramatic reduction in the number of operations required to perform the update of the kinetic equation. In doing so, it may be possible to reduce the scaling of the algorithm from to for the configuration space component of the update and to for the velocity space component of the update. In this case, is the number of basis functions in the expansion in configuration space for the DG discretization of Maxwell’s equations. In fact, we find in a preliminary test of a modal basis expansion of the Vlasov-Maxwell system that these scaling do hold, encouraging further study of this alternative formulation. Traditionally, this approach is not justified because of how expensive the modal surface updates can become in 2D and 3D, but our preliminary tests additionally reveal a turnover from 3D to 4D, wherein even the surface integrals are much cheaper with the orthogonal modal basis. Other improvements to the solver can come from hardware considerations. Current computing paradigms are quickly moving towards many-core devices, which dramatically increase the number of floating point operations the architecture can do per second at the cost of available memory. This paradigm favors arithmetically intense algorithms like DG, as long as one embraces hybrid parallel programming models so that memory can be shared amongst cores and thus shared amongst floating point units.
Finally, we consider further additions to the solver beyond the core algorithm which would further improve the tool’s capabilities. As mentioned in the introduction, the algorithm derived and implemented here is in both the non-relativistic and collisionless limit of the Vlasov-Maxwell system of equations. Of course, there exist many plasma systems, especially in astrophysics, where relativistic effects are important, such as gamma ray bursts or the magnetospheres of neutron stars. And many plasma systems, such as the solar wind, are not in fact collisionless, but weakly collisional. The inclusion of the Lorentz boost factors and a discretized collision operator, either reduced such as the Bhatnagar-Gross-Krook[32], Eq. (7), or Lenard-Bernstein, Eq. (8), operators [33, 34, 35], or the full Landau-Fokker-Planck, Eq. (4), operator [31], would be useful extensions to the model. It should be noted that the inclusion of relativistic effects in the Vlasov equation necessitates the inclusion of relativistic effects in the collision operator, e.g., the implementation of the Beliav-Budker operator [61] instead of the Landau-Fokker-Planck operator, or some relativistic generalization of the example reduced collision operators. While the inclusion of the physics of relativistic effects would be vital for simulating extreme astrophysical environments such as pulsar wind nebulae, just as the study of dissipation and entropy production requires a discrete collision operator, there are challenges to the inclusion of these model extensions which are worth acknowledging. For example, when adding relativistic effects, one key element to our algorithm that must be respected is that the integrals in the discrete-weak form must be evaluated exactly. Unlike other systems in which DG has found success, such as the solutions of fluid equations, where aliasing errors from under-integration are often tolerated, no such concession can be made for a DG discretization of the Vlasov equation, lest one is willing to risk numerical instabilities from a lack of energy and particle conservation. Unfortunately, as can be seen in Eq. (2) and Eq. (3), the Lorentz boost factors in the relativistic Vlasov equation will inevitably lead to aliasing errors in our current algorithm due to the fact that standard Gauss-Legendre quadrature will fail to exactly evaluate the square roots and rational functions which arise from a DG discretization of the relativistic Vlasov equation. Even an algorithmic change to pre-evaluate the integrals, such as the one we are pursuing with the modal basis, may not work as closed form expressions for the integrals may not exist. The inclusion of a discrete collision operator has its own challenges, namely that collisions tend to manifest as a stiff source term in our equations, similar to the issues seen in the solution of advection-diffusion equations [62, 63]. This challenge would normally motivate an operator split so as to treat the collision term implicitly, but since collision operators such as the Lenard-Bernstein and Landau-Fokker-Planck operators are nonlinear, such a split is nontrivial and may require very sophisticated implicit solves. These issues, while challenging, do have solutions in specific, non-general cases. For example, a worthwhile avenue to pursue for including relativistic effects would be to Taylor expand the Lorentz Boost factors locally in a momentum cell, thus allowing the treatment of some relativistic effects while maintaining an alias-free algorithm. And since many plasma systems are weakly collisional, a low enough collisionality such that the collision operator is not as stiff as the advection term could be justified for certain applications. These solutions would very likely break down when attempting to model highly relativistic, highly collisional plasmas; however, a path forward exists for including the physics of collisions and special relativity without solving the fully general problem.
Beyond the addition of a physical model for collisions and the inclusion of relativistic effects, there is much to be gained in added flexibility of our velocity space discretization. Currently, velocity space boundary conditions remain an unphysical characteristic of the discrete system. An assumption that no particles will reach velocities at the edge of velocity space is a dangerous one, as sufficient particle acceleration or heating can quickly lead to numerical issues when the plasma interacts with the velocity space “wall.” Extending velocity space to via something like a “super-cell” wherein the volume term is integrated with Gauss-Laguerre instead of Gauss-Legendre quadrature, as some gyrokinetic codes currently employ [64, 65], would allow for the treatment of plasma systems which experience non-trivial amounts of particle acceleration and/or heating. Non-uniform grids in velocity space would allow one to concentrate velocity resolution around expected resonance processes. Adaptive mesh refinement in velocity space [66, 67] would prove invaluable for highly nonlinear simulations such as the turbulence calculation presented in this paper, as velocity space requirements vary dramatically throughout the domain, from minimal in quiescent parts of the plasma to quite demanding in the middle of reconnecting current sheets such as the one which forms in the center of the domain in the Orszag-Tang vortex. Ideally, such mesh refinement would be anisotropic, as a regridding procedure in five or six dimensions would be quite demanding, and configuration space mesh refinement is not strictly necessary; for many problems of interest, the configuration space scales one wants to resolve are well-defined. Such flexibility in velocity space grids is not guaranteed to provide significant advantage over uniform, fixed, grids, as processes such as phase-mixing can create structure everywhere in velocity space at arbitrarily small velocity space scales; however, for simulations which include significant particle acceleration such as the electrostatic shock in Section 5.2.3, the generation of the highly energetic tail in the proton distribution function would be much cheaper. And relativistic laser-plasma Vlasov codes in 1X1P have demonstrated significant gains with AMR [15]. Most importantly, as discussed in Sections 5.2.1 and 5.3.2, there is currently nothing inhibiting the distribution function from going negative, nor is there anything preventing the equation system from accumulating errors in the divergence relations in Maxwell’s equations, Eqns. 11 and 12. Although the fluid moments such as density and temperature are positive definite for the extent of the run, the distribution function likewise should also be positive definite, and although the divergence errors remain small, charge should be conserved and the divergence of the magnetic field should remain zero. The implementation of a positivity, and perhaps monotonicity, enforcing limiter, either to the solution or the differential operator, which preserves particles and energy, and the implementation of a curl-conforming DG method for Maxwell’s equations, such as the one outlined in Pintos et al. [60], would be the two most valuable improvements of all the algorithmic enhancements discussed here.
Nonetheless, despite many possibilities for future development, we have demonstrated the robustness of a DG Vlasov-Maxwell solver through a variety of benchmarks. The level of signal present in the continuum discretization of the Vlasov-Maxwell system is quite substantial. The Orszag-Tang run exhibits greater than eight orders of magnitude of signal in the distribution function. Given cutting edge diagnostics such as the field particle correlations [48, 49, 50], the value of having a noise-free distribution function is made manifest. Many exciting problems can be tackled with this tool, spanning the space of plasma systems found in the lab, such as radio-frequency heating problems, and space and astrophysical plasmas, such as ongoing attempts to quantify dissipation in plasma turbulence in the heliosphere.
Acknowledgements
The authors acknowledge fruitful discussions with Greg Hammett, Kristopher Klein, and Petr Cagas. The work of J. Juno and J. TenBarge was supported by the National Science Foundation SHINE award No. AGS-1622306 and the UMD DOE grant DE-FG02-93ER54197. The work of W. Dorland was also supported by the UMD DOE grant DE-FG02-93ER54197. The work of A. Hakim and E. Shi was supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466, through the Max-Planck/Princeton Center for Plasma Physics and the Princeton Plasma Physics Laboratory. A. Hakim was also supported by the Air Force Office of Scientific Research under Grant No. FA9550-15-1-0193. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.
References
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] C. Birdsall, A. B. Langdon, Plasma Physics Via Computer Simulation, Institute of Physics Publishing, 1990.
- 2[2] E. A. Frieman, L. Chen, Nonlinear gyrokinetic equations for low-frequency electromagnetic waves in general plasma equilibria, Phys. Fluids 25 (1982) 502–508.
- 3[3] A. J. Brizard, T. S. Hahm, Foundations of nonlinear gyrokinetic theory, Rev. Mod. Phys. 79 (2007) 421–468. doi:10.1103/Rev Mod Phys.79.421 . · doi ↗
- 4[4] E. Camporeale, G. L. Delzanno, B. K. Bergen, J. D. Moulton, On the velocity space discretization for the Vlasov-Poisson system: Comparison between implicit Hermite spectral and Particle-in-Cell methods, Comp. Phys. Comm. 198 (2016) 47–58.
- 5[5] M. W. Kunz, J. M. Stone, X.-N. Bai, Pegasus: A new hybrid-kinetic particle-in-cell code for astrophysical plasma dynamics, J. Comp. Phys. 259 (2014) 154–174. doi:10.1016/j.jcp.2013.11.035 . · doi ↗
- 6[6] F. Valentini, P. Trávníček, F. Califano, P. Hellinger, A. Mangeney, A hybrid-Vlasov model based on the current advance method for the simulation of collisionless magnetized plasma, J. Comp. Phys. 225 (2007) 753–770. doi:10.1016/j.jcp.2007.01.001 . · doi ↗
- 7[7] F. Valentini, A. Iazzolino, P. Veltri, Numerical study of ion-cyclotron resonant interaction via hybrid-Vlasov simulations, Phys. Plasmas 17 (5). doi:10.1063/1.3420278 . · doi ↗
- 8[8] A. Greco, F. Valentini, S. Servidio, W. H. Matthaeus, Inhomogeneous kinetic effects related to intermittent magnetic discontinuities, Phys. Rev. E 86 (6) (2012) 1–6. doi:10.1103/Phys Rev E.86.066405 . · doi ↗