Sparse Grid Discontinuous Galerkin Methods for the VlasovMaxwell System
Abstract
In this paper, we develop sparse grid discontinuous Galerkin (DG) schemes for the VlasovMaxwell (VM) equations. The VM system is a fundamental kinetic model in plasma physics, and its numerical computations are quite demanding, due to its intrinsic highdimensionality and the need to retain many properties of the physical solutions. To break the curse of dimensionality, we consider the sparse grid DG methods that were recently developed in [20, 21] for transport equations. Such methods are based on multiwavelets on tensorized nested grids and can significantly reduce the numbers of degrees of freedom. We formulate two versions of the schemes: sparse grid DG and adaptive sparse grid DG methods for the VM system. Their key properties and implementation details are discussed. Accuracy and robustness are demonstrated by numerical tests, with emphasis on comparison of the performance of the two methods, as well as with their full grid counterparts.
100
Keywords: discontinuous Galerkin methods; sparse grids; VlasovMaxwell system; streaming Weibel instability; Landau damping.
1 Introduction
The VlasovMaxwell (VM) system is a fundamental model in plasma physics for describing the dynamics of collisionless magnetized plasmas, which finds diverse applications in science and engineering, including thermonuclear fusion, satellite amplifiers, highpower microwave generation, to name a few. In this paper, we study the VM system that describes the evolution of a single species of nonrelativistic electrons under the selfconsistent electromagnetic field while the ions are treated as uniform fixed background. Under the scaling of the characteristic time by the inverse of the plasma frequency , length by the Debye length , and electric and magnetic fields by (with the electron mass, the speed of light, and the electron charge), the dimensionless form of the VM system is
(1a)  
(1b)  
(1c)  
(1d) 
with
where the equations are defined on denotes position in physical space, and which is the velocity space. Here is the distribution function of electrons at position with velocity at time , is the electric field, is the magnetic field, is the electron charge density, and is the current density. The charge density of background ions is denoted by , which is chosen to satisfy total charge neutrality, . Ideally however numerical computation requires a truncation of the space and the assumption that is compactly supported on In this paper, for simplicity, we will assume to be boxshaped domains.
The simulations of VM systems are quite challenging. Particleincell (PIC) methods [5, 22] have long been very popular numerical tools, mainly because they can generate reasonable results with relatively low computational cost. However, as a MonteCarlo type approach, the PIC methods are known to suffer from the statistical noise, which is with being the number of sampling particles. Such an inherent low order error of PIC methods prevents accurate description of physics of interest, when, for instance, the tail of the distribution function needs to be resolved accurately. In recent years, there has been growing interest in deterministic simulations of the Vlasov equation. In the deterministic framework, the schemes are free of statistical noise and hence able to generate highly accurate results in phase space. In the context of VM simulations, Califano et al. have used a semiLagrangian approach to compute the streaming Weibel (SW) instability [10], current filamentation instability [28], magnetic vortices [9], magnetic reconnection [8]. Also, various methods have been proposed for the relativistic VM system [33, 4, 34, 23]. In this paper, we are interested in a class of successful deterministic Vlasov solvers based on the discontinuous Galerkin (DG) finite element discretization, because of not only their provable convergence and accommodation for adaptivity and parallel implementations, but also their excellent conservation property and superior performance in long time wavelike simulations. Those distinguishing properties of DG methods are very much desired for the VM simulations, and they have been previously employed to solve VM system [13, 12] and the relativistic VM system [37]. However, due to the curse of dimensionality, traditional deterministic approaches including the DG methods are not competitive for high dimensional Vlasov simulations, even with the aid of high performance computing systems.
To break the curse of dimensionality, this paper will focus on the sparse grid approach. The sparse grid method [7, 17] has long been an effective numerical tool to reduce the degrees of freedom for highdimensional grid based methods. In the context of wavelets or sparse grid methods for kinetic transport equations, we mention the work of using waveletMRA methods for Vlasov equations [3], the combination technique for linear gyrokinetics [27], sparse adaptive finite element method [36], sparse discrete ordinates method [18] and sparse tensor spherical harmonics [19] for radiative transfer, among many others. In [35, 20], a class of sparse grid DG schemes were proposed for solving highdimensional partial differential equations (PDEs) based on a novel sparse DG finite element approximation space. Such a sparse grid space can be regarded as a proper truncation of the standard tenor approximation space, which reduces degrees of freedom of from to , where is the uniform mesh size in each direction and is the dimension of the problem. Coupling the DG weak formulation with the sparse grid space, the resulting sparse grid DG method is demonstrated to save computational and storage cost without deteriorating the approximation quality by much. In particular, when applied to a dimensional linear transport equation, the scheme is proven to be stable and convergent with rate in norm for a smooth enough solution, where is the degree of polynomials [20]. Motivated by the development of sparse grid DG method [20] and the adaptive multiresolution DG method [21] for transport equations, it is of interest to this paper to develop sparse grid DG methods for solving the VM system. The proposed methods are well suited for VM simulations, due to their ability to handle high dimensional convection dominated problems, the ability to capture the main structures of the solution with feasible computational resource and the overall good performance in conservation of physical quantities in long time simulations.
2 Numerical methods
In this section, we describe two sparse grid DG methods for the VM system: the standard sparse grid DG method and the adaptive sparse grid DG method. We first review the finite element space on sparse grid introduced in [35], and then describe the details of the schemes when applied to the VM system.
2.1 DG finite element space on sparse grid
In this subsection, we review the notations of DG finite element space on sparse grid. First, we introduce the hierarchical decomposition of piecewise polynomial space in one dimension. Without loss of generality, consider the interval , we define a set of nested grids, where the th level grid consists of uniform cells , for any The nested grids result in the nested piecewise polynomial spaces. In particular, let
be the usual piecewise polynomials of degree at most on the th level grid . Then, we have
We can now define the multiwavelet subspace , as the orthogonal complement of in with respect to the inner product on , i.e.,
Note that the space , for , represents the finer level details when the grid is refined and this is the key to the reduction of the degrees of freedom in higher dimensions. We further let , which is standard piecewise polynomial space of degree on . Therefore, we have found a hierarchical representation of the standard piecewise polynomial space on as . In [35], we used the orthonormal basis functions for space which are constructed based on the onedimensional orthonormal multiwavelet bases first introduced in [1]. We refer readers to [1] and [35] for more details. We now denote the basis functions in level as
and they are orthonormal, i.e.
Below we give an example of the basis functions. We denote the bases for as
where are functions supported on depending on . For instance, when , the formulas for restricted to the interval are
The functions are extended to as even or odd functions according to the parity of :
and are zero outside . Then, the bases for are defined as
Now we are ready to prescribe the sparse finite element space in dimensions on . Notice similar discussions apply for any boxshaped domain in dimensions. First we recall some basic notations about multiindices. For a multiindex , where denotes the set of nonnegative integers, the and norms are defined as
The componentwise arithmetic operations and relational operations are defined as
By making use of the multiindex notation, we indicate by the mesh level in a multivariate sense, where denotes nonnegative integers. We consider the tensorproduct rectangular grid with mesh size Based on the grid , an elementary cell is denoted by , and
is the tensorproduct piecewise polynomial space, where denotes polynomials of degree up to in each dimension on cell . If , the grid and space will be denoted by and , respectively. For the increment space the orthonormal basis functions can be defined as
where denote orthonormal bases in mth dimension defined in onedimensional case. With the onedimensional hierarchical decomposition, we have
The sparse finite element approximation space on mesh we use in this paper, on the other hand, is defined by [35, 20]
This is a subset of , and its number of degrees of freedom scales as or , where denotes the finest mesh size in each direction [35]. This is significantly less than that of with number of elements when is large. The sparse grid DG scheme in Section 2.2 is defined using this space. The adaptive sparse grid DG scheme, however, relies on an adaptive choice of the space, and will be discussed in details in Section 2.3.
2.2 The sparse grid DG scheme
Below, we formulate a DG scheme in the sparse finite element space for the VM system (1a)(1b) inspired by [13]. On the PDE level, the two equations in (1c) involving the divergence of the magnetic and electric fields can be derived from the remaining part of the VM system. However, how to impose (1c) numerically is an important and nontrivial issue [29, 2, 24]. This will be studied in our future work.
Using the notations introduced in the previous subsections, and let and be the dimension of and , respectively, we consider the partitions of the domain into mesh level in all directions. We distinguish between the  and directions. Let and be the collection of all elements in and , respectively. Let be the union of the edges for all elements in , similarly let be the union of the edges for all elements in . Furthermore, with and being the union of interior and boundary edges of , respectively.
For piecewise functions, we further introduce the jumps and averages as follows. Suppose and are two elements in . For any edge , with as the outward unit normal to , , and , the jumps across are defined as
and the averages are
When considering periodic boundary conditions, the jumps and averages can be naturally defined on the boundary edges.
By replacing the subscript with , the jumps and averages can be defined similarly for the direction. For a boundary edge with being the outward unit normal, we use
(2) 
This is consistent with the fact that the exact solution is assumed to be compactly supported in domain.
We are now ready to describe the scheme. The sparse discrete spaces on and we use are defined as
Following [20], the semidiscrete DG methods for the VM system are: to find , , such that for any , ,
(3a)  
(3b)  
(3c) 
with
(4) 
All “hat” functions are numerical fluxes. For the Vlasov part, we adopt the global LaxFriedrichs flux:
(5a)  
(5b) 
where and , where the maximum is taken for all possible and at time in the computational domain. For the Maxwell part, we use the upwind flux
(6a)  
(6b) 
or the alternating fluxes
(7) 
Previous studies in [12, 13] have demonstrated the importance of numerical flux on the quality of DG simulations. It was understood that a dissipative flux choice is desired for the Vlasov equation. The alternating and upwind fluxes are both optimal in order for the Maxwell solver, but alternating flux is energyconserving for the Maxwell’s equation, while upwind flux is not. We will not consider the central flux based on the recommendations made in [13].
Below, we will summarize the conservation and stability properties of the semidiscrete scheme. The proof is very similar to those in [13] and is thus omitted.
Theorem 2.1 (Mass conservation)
The numerical solution with satisfies
(8) 
where
Equivalently, with , for any , the following holds:
(9) 
Theorem 2.2 (Energy conservation)
In the theorems above, terms like come from numerical errors of the velocity boundary truncation from an infinite to a finite domain. If those terms are negligible, i.e. when is chosen sufficiently large, we can conclude that the scheme is massconservative, and energyconservative if the alternating flux is used for the Maxwell solver.
Theorem 2.3 (stability of )
For , the numerical solution satisfies
As for time, we use the total variation diminishing RungeKutta (TVDRK) methods [32] to solve the ordinary differential equations resulting from the discretization (3a)(3c), denoted by In particular, we use the following thirdorder TVDRK method in this paper
(10)  
where represents the numerical solutions at time level .
Finally, we mention some details about implementation. A key to computational efficiency is the efficient evaluation of multidimensional integrations. To compute multidimensional integration, we apply the unidirectional principle. For example, if we want to evaluate with , it is equivalent to multiplication of onedimensional integrals Based on the hierarchical structure of the basis functions, we only need some small overhead to compute onedimensional integrals and assemble them to obtain the multidimensional integrations. In (3a)(3c), the numerical integrations with coefficients and (which belong to ) can be done very efficiently because they can all be computed using this trick. This procedure is performed one time before time evolution starts, and the sparsity of the matrices due to the multiwavelet basis structures is utilized to accelerate the computation [35].
2.3 The adaptive sparse grid DG scheme
In this subsection, we will describe the adaptive sparse grid DG (also called adaptive multiresolution DG) method [21] for the VM system. The motivation to study the adaptive scheme is to offer better numerical resolutions for the probability density function. It is known that the solution to Vlasov equation develops filamentation, therefore the standard sparse grid method can not offer enough resolution when becomes large (see [21] for comparison in the case of the VlasovPoisson system). Therefore, in this paper, we also consider the adaptive scheme.
The main idea of the algorithm is not to use in a predetermined fashion, but rather to choose a subspace of adaptively. In this work, adaptivity is only implemented for not for the lowerdimensional variables which are computed using the full finite element space. This will not cause much additional cost because the Maxwell’s equations are in lower dimension than the Vlasov equation. If the computational and storage cost is a big concern, then the adaptive strategy can be potentially applied to the Maxwell’s part as well.
The algorithm is described in details below. Given a maximum mesh level and an accuracy threshold , we first use the adaptive multiresolution projection algorithm [21] to get the numerical initial condition for the DG scheme. For completeness of the paper, the details of the method is listed below in Algorithm 1.
Algorithm 1: Adaptive multiresolution projection
Input: Function .
Parameters: Maximum level polynomial degree error threshold
Output: Hash table , leaf table and projected solution

Project onto the coarsest level of mesh, e.g., level 0. Add all elements to the hash table (active list). Define an element without children as a leaf element, and add all the leaf elements to the leaf table (a smaller hash table).

For each leaf element in the leaf table, if the refinement criteria
(11) holds, then we consider its child elements: for a child element , if it has not been added to the table , then compute the detail coefficients and add to both table and table . For its parent elements in , we increase the number of children by one.

Remove the parent elements from table for all the newly added elements.

Repeat step 2  step 3, until no element can be further added.
When the adaptive projection algorithm completes, it will generate a hash table , leaf table and . We denote the approximation space and it is a subspace of On the other hand, we compute the initializations of by a simple projection of onto .
Then we begin the time evolution algorithm which consists of several key steps. The first step is the prediction step, which means that, given the hash table that stores the numerical solution at time step and the associated leaf table , we need to predict the location where the details becomes significant at the next time step , then add more elements in order to capture the fine structures. We solve for from to based on the solutions and at time . The forward Euler discretization is used as the time integrator in this step and we denote the predicted solution at by
The second step is the refinement step according to . We traverse the hash table and if an element satisfies the refinement criteria (11), indicating that such an element becomes significant at the next time step, then we need to add its children elements to and provided they are not added yet, and set the associated detail coefficients to zero. We also need to make sure that all the parent elements of the newly added element are in (i.e., no “hole” is allowed in the hash table) and increase the number of children for all its parent elements by one. This step generates the updated hash table and leaf table . Note that (11) corresponds to the norm based refinement criteria in [21].
Based on the updated hash table , the third step is to evolve the numerical solution by the DG weak formulation with space . Namely, we solve for from to , by the TVDRK scheme (10) to generate the precoarsened numerical solution . Meanwhile, we also evolve the numerical solutions and from to with space .
The last step is to coarsen the mesh by removing elements that become insignificant at time level The hash table that stores the numerical solution is recursively coarsened by the following procedure. The leaf table is traversed, and if an element satisfies the coarsening criterion
(12) 
where is a prescribed error constant, then we remove the element from both table and , and set the associated coefficients . For each of its parent elements in table , we decrease the number of children by one. If the number becomes zero, i.e, the element has no child any more, then it is added to the leaf table accordingly. Repeat the coarsening procedure until no element can be removed from the table . The output of this coarsening procedure are the updated hash table and leaf table, denoted by and respectively, and the compressed numerical solution . In practice, is chosen to be smaller than for safety. In the simulations presented in this paper, we use . For the adaptive sparse grid DG scheme, the properties of conservation of moments are not as clear as the sparse grid DG schemes, and will be subject to the error thresholds , see [21].
3 Numerical results
In this section, we consider two numerical tests to benchmark the performance of the proposed schemes.
3.1 1D2V streaming Weibel instability
Here, we consider the streaming Weibel (SW) instability, which has been previously considered both analytically and numerically [10, 13, 12].
Example 3.1
In this case, a reduced version of the single species VM system with one spatial variable , and two velocity variables is considered:
(13)  
(14) 
where
(15) 
Here, is the distribution function of electrons, is the 2D electric field, and is the 1D magnetic field. Ions are assumed to form a constant background.
The initial condition is given by
(16)  
(17) 
where denotes the amplitude of the initial perturbation to the magnetic field, is the thermal velocity, which is take to be , and is a parameter measuring the symmetry of the electron beams. Note that when , this initial condition is an equilibrium state representing counterstreaming beams propagating perpendicular to the direction of inhomogeneity. As in [10], the instability can be triggered by taking as a perturbation of the initial magnetic field. The computational domain is chosen to be , where and . We consider two sets of parameters as in [10]
which lead to initially symmetric and strongly nonsymmetric counterstreaming electron beams, respectively. For all numerical simulations, unless otherwise noted, the time step is chosen according to the mesh on the most refined level, i.e.
where is the maximum wave propagation speed in th direction, and we take .
Accuracy test and comparisons. It is wellknown that the VM system is time reversible, and such property provides practical means to test accuracy for VM solvers. To elaborate, let be the initial condition and be the solution at for the VM system. When we reverse the velocity field of the solution and the magnetic filed, yielding , and evolve the VM system again to , then we can recover , which is the initial condition with the reverse velocity field and the reverse magnetic field. For the accuracy test, we compute the solutions to and then reversely back to , and compare the numerical solutions with the initial condition. The errors of are defined as
for :  
for :  (18)  
for : 
where are the numerical approximations to the exact solutions .
We test accuracy for the sparse grid DG method with on different levels of meshes. For , we take to match the temporal and spatial orders in the convergence study. The errors and orders of the numerical solutions with upwind and alternating fluxes for parameter choice 1 are reported in Tables 1 and 2. Similar to [20], we observe at least th order accuracy for both fluxes. When comparing the results of the two tables, we find that the errors are identical for , while there are some slight differences in the errors of the electromagnetic fields. This is expected because the solvers only differ in the discretization of Maxwell’s equations.
error  order  error  order  error  order  

7  1.25E01  7.49E08  9.82E08  
8  4.44E02  1.49  2.86E08  1.39  2.32E08  2.08  
9  2.01E02  1.14  5.01E09  2.51  5.61E09  2.05  
10  5.74E03  1.81  1.50E09  1.74  1.32E09  2.09  
7  3.09E02  1.37E09  9.25E10  
8  6.65E03  2.22  1.40E10  3.29  9.09E11  3.35  
9  2.56E03  1.39  2.57E11  2.45  9.13E12  3.32  
10  4.07E04  2.65  4.67E12  2.46  1.24E12  2.88  
7  8.47E03  1.45E11  7.82E12  
8  9.37E04  3.18  1.18E12  3.62  3.48E13  4.49  
9  9.48E05  3.31  7.09E14  4.06  3.87E14  3.17  
10  8.48E06  3.48  9.49E15  2.90  1.40E15  4.79 
error  order  error  order  error  order  

7  1.25E01  7.67E08  6.32E08  
8  4.44E02  1.49  2.96E08  1.37  1.11E08  2.51  
9  2.01E02  1.14  5.83E09  2.34  2.36E09  2.23  
10  5.74E03  1.81  1.83E09  1.67  3.53E10  2.74  