High Order and Energy preserving Discontinuous Galerkin Methods for the VlasovPoisson system
Abstract.
We present a computational study for a family of discontinuous Galerkin methods for the one dimensional VlasovPoisson system, recently introduced in [acs0]. We introduce a slight modification of the methods to allow for feasible computations while preserving the properties of the original methods. We study numerically the verification of the theoretical and convergence analysis, discussing also the conservation properties of the schemes. The methods are validated through their application to some of the benchmarks in the simulation of plasma physics.
Key words and phrases:
plasma physics, discontinuous Galerkin, Vlasov Poisson system, energy conservation1991 Mathematics Subject Classification:
82C80, 65M60, 65M12, 82A70Numerical simulation has become a major tool for understanding the complex behavior of a plasma or a particle beam in many situations. This is due not only to the large number of physical applications and technological implications of the behavior of plasmas, but also to the intrinsic difficulties of the models used to describe such behavior. In fact, it was recognized long time ago that there does not exist any fully satisfactory macroscopic model (fluid equations) which can be used to describe the particle interaction in laserfusion problems. In contrast, microscopic models (kinetic equations) can provide a more accurate description of the plasmas.
One of the simplest model problems that is currently used in the simulation of plasmas is the VlasovPoisson system. Such system describes the evolution of a plasma of charged particles (electrons and ions) under the effects of the transport and selfconsistent electric field. The unknown, typically denoted by f(x,v,t) (with x standing for position, v for velocity and t for time), represents the distribution function of particles (ions, electrons, etc.) in the phase space. The coupling with a selfconsistent electrostatic field (neglecting magnetic effects) is taken into account through the Poisson equation. The nonlinear structure of the system prevents from obtaining analytical solutions, except for a few academic cases (see the surveys [glassey, Bouchut, dolbeault] for a good description on the state of the art of the mathematical analysis of the problem). Therefore, numerical simulations have to be performed to study realistic physical phenomena.
At the present time, there can be distinguished two main classes of numerical methods for simulating plasmas; Lagrangian (or probabilistic) and Eulerian (or deterministic) methods. The former class include all different types of particle methods [BirdLang, CoRa84, WONumer96, GanVic89, pic1, pic3, pic2] and has been a preferred choice since the beginnings of numerical simulations in plasma physics in the 60^{\prime}s, due to their simplicity and low computational cost. The basic idea behind these methods is to approximate the motion of the plasma by a finite number of macroparticles in the phase space whose trajectories are computed from the characteristics of the Vlasov equation, while the electrostatic field is computed by collecting the charge density on a fixed mesh of the physical space. Although this class of methods represents a feasible option and potentially might allow for resolving the whole 3+3+1 dimensional problem, their inherent numerical noise precludes from obtaining an accurate description of the distribution function in
the phase space in many interesting cases. This lack of precision can be overcome by using a method from the second class; an Eulerian solver. These type of methods are nothing but classical (or new) numerical schemes discretizing the Vlasov equation on a (fixed) mesh of the phase space. Among them, the most widely used are finite volumes [Fijal99, filbet] and semilagrangian methods [filbet2, Bes04, mcp0, CV07, crson00, sl1, sl2, sl3]. Finite volumes (FV) are a simple and inexpensive option, but in general are low order if one wants to retain the basic conservative properties of the scheme.
Semilagrangian schemes (sometimes consider inbetween Eulerian and Lagrangian solvers) have become a popular option, since they can achieve high order allowing at the same time for time integration with large time steps. However, special care in needed to compute the origin of the characteristics with high order interpolation without spoiling the local character of the reconstruction. A nice numerical study comparing some of the different methods use in plasma simulations is presented in [comparison].
In this paper we present a computational study with discontinuous Galerkin (DG) methods. DG methods are finite element methods that use discontinuous polynomials. Their local construction endow the methods with good local conservation properties without sacrificing the order of accuracy. This is one of the main motivations for their use in plasma simulations. But it also provides the methods a builtin parallelism which allows for parallelization of the resulting algorithms. The methods have also many other attractive features: they are extremely flexible in handling hpadaptivity, the boundary conditions are imposed weakly and the DG mass matrices are blockdiagonal which results in very efficient timestepping algorithms in the context of timedependent problems, as it is the case here. In spite that nowdays, DG methods are consider for approximating problems of almost any kind, their use for kinetic equations, and more particularly for simulation of plasmas has only been contemplated very recently. Among the computational works, we mention the use of DG in a multiwaterbag approximation of the VP system in [BBBB09]; a piecewise constant DG solver for VP in [gamba0] (which require extremely fine meshes) and semilagrangian schemes combined with high order DG interpolation are presented in [seal, qiuCW]. In both works, the authors also use the positivity preserving limiter introduced in [zsposit0].
A theoretical work has been presented in [acs0, acs1], where the authors have introduced and analyzed a family of semidiscrete DG schemes for the VP system with periodic boundary conditions, for the one and multidimensional cases, respectively. The authors show optimal error estimates for both the distribution function and the electrostatic field, and they study the conservation properties of the proposed schemes. Due to the local construction of the DG schemes, total mass (or charge) conservation is shown to hold easily. This property is essential in the numerical approximation to VP, since it is required for guaranteeing the wellposedness of the related Poisson problem. The authors also propose a novel DG scheme that is shown to preserve the total energy of the VP system. Their proof however requires the assumption that the DG finite element spaces contain at least all quadratic polynomials.
In this work, we undertake the issues of verification and validation of these family of DG schemes, for the onedimensional VP system. To accomplish both tasks, we first discuss how the schemes can be efficiently implemented in practice, even in parallel. For the space discretization, the methods introduced in [acs0] are based on the coupling of a DG discretization for the Vlasov equation (transport equation) together with a mixed finite element (possibly discontinuous) approximation to the Poisson problem. Here, however, we present two slight variations of the DG approximation for the Vlasov equation, to allow for feasible computations. The modifications are done in the definition of the numerical flux involving the coupling with the approximate electrostatic field (hence the nonlinearity). The definition in [acs0] would require the computation, at each time step of the zeros of the approximate electrostatic field, which would increase substantially the cost, taking into account that we use high order approximations. Nevertheless, as we show here, the slight variation in the schemes does not affect the optimal accuracy of the methods. Furthermore, since the new definition of the flux is still consistent, the mass and energy conservation can still be guaranteed (even at the theoretical level). Also, here we demonstrate numerically that for the energy preserving scheme given in [acs0], it is indeed necessary (and not a technical restriction due to the proof) to use finite element spaces spaces containing all quadratic polynomials.
For the time discretization we stick to a simple fourth order explicit Runge Kutta (RK) method, the socalled RK4 or classic RungeKutta [hairer0]. The reason for not using total variation diminishing (TVD) RK integrator is twofold. On the one hand, in our simulations we have observed no numerical evidence of any essential benefit of the TVD integrator over the standard RK method (probably due to the smoothness of the solution). On the other hand, since we focus on high order methods (for the space discretization), the time integration should be accomplished also with some high order time integration scheme. As is well known [tvd0], a fourth (or higher) order TVD RK, would require for the computation of the internal stages, the evaluation of the operator and its adjoint, due to the presence of some negative coefficients in the corresponding TVDRK tableau. This would substantially increase the cost (and storage) of the overall procedure, without any significant benefit. With the fully discrete schemes, we verify numerically the theory developed in [acs0]; both the error analysis together with the conservation properties.
The second goal of the paper is to validate the methods by studying their performance in approximating some of the classical benchmark problems in plasma physics. Here we consider the linear and nonlinear Landaudamping together with two benchmarks related to the two stream instability problems. We compare our numerical results with those available in literature, getting always at least the same outcomes. In particular, we show the benefit of using the energy preserving high order DG method for the numerical simulations (since no extremely refined meshes are needed and the code can be parallelized).
In the last part of the paper, we consider the application of the schemes for the boundary value problem of the VP system studied in [filbet_shu].
This problem models the evolution of a collisionless electron gas under the influence of a electrostatic field E in an interval [0,1], with electrons emitted at one end and absorbed at the other end of the interval. Due to the absorbing boundary condition, it has been proved theoretically the distribution function f might become discontinuous in finite time, depending on the sign and magnitude of the electrostatic field at the boundary.
Although the DG methods we consider in this paper were not originally designed to approximate such problem, we study here the ability of the methods to capture the discontinuity.
The results however, are not completely satisfactory, since we do not always (at all the times) capture the behaviour of the solution predicted by the theory developed in [filbet_shu]. A possible reason is the weak nature of the singularity, but it might also happen that as the time evolves the full discretization is adding too much artificial viscosity, which does not allow the methods to capture completely the singularity. This issue together with the tuning of the schemes to capture correctly the singularities (at all times) will be the subject of future study.
The outline of the paper is as follows. In Section 1 we describe the main properties of the continuous problem and introduce the basic notations related to the discrete DG methods. We then introduce the numerical methods we consider discussing also their main properties in Section 2. In Section LABEL:sec:4 we consider the full discretization and deal with the implementation issues related to the schemes. Section LABEL:sec:6 is devoted to the validation and convergence study of the schemes. We present extensive numerical tests and consider the application of the methods for the simulation of Landau damping and two different tests related to the nonlinear two stream instability. In section LABEL:sec:7 we examine the application of the considered DG methods for approximating a VlasovPoisson boundary value problem (no periodic boundary conditions). Finally, we derive some conclusion in section LABEL:sec:fin
Notation: Throughout this paper, we use the standard notation for Sobolev spaces (see [Adams75]). For a bounded domain B\subset\mathbb{R}^{2}, we denote by H^{m}(B) the standard Sobolev space of order m\geq 0 and by \\cdot\_{m,B} and \cdot_{m,B} the usual Sobolev norm and seminorm, respectively. For m=0, we write L^{2}(B) instead of H^{0}(B). We shall denote by H^{m}(\mathcal{I})/\mathbb{R} the quotient space consisting of equivalence classes of elements of H^{m}(\mathcal{I}) differing by constants; for m=0 it is denoted by L^{2}(\mathcal{I})/\mathbb{R}. We shall indicate by L^{2}_{0}(B) the space of L^{2}(B) functions having zero average over B.
1. The VlasovPoisson system and basic notation
In this section we introduce the Vlasov Poisson system and recall some of its properties. In the last part of the section, we also introduce the basic notation required for describing the numerical methods we consider.
1.1. Continuous problem: the VlasovPoisson system
We consider a noncollisional plasma of charged particles (electrons and ions). For simplicity, we assume that the properties of the plasma are one dimensional and we take into account only the electrostatic forces, thus neglecting the electromagnetic effects. We denote by f=f(x,v,t) the electron distribution function and by E(x,t)=\Phi_{x}(x,t) the electrostatic field. The VlasovPoisson equations of the plasma in dimensionless variables can be rewritten as,
(1.1)  \displaystyle\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}\Phi% _{x}\frac{\partial f}{\partial v}  \displaystyle=0  \displaystyle(x,v,t)\,\,\in\,\,\Omega_{x}\times\mathbb{R}\times[0,t_{f}],  
(1.2)  \displaystyle\Phi_{xx}  \displaystyle=\rho(x,t)1  \displaystyle(x,t)\,\,\in\,\,\Omega_{x}\times[0,t_{f}], 
where v denotes the velocity of the charged particles and \rho(x,t) is the charge density defined by
\rho(x,t)=\displaystyle{\int_{\mathbb{R}}f(x,v,t)dv}\quad\forall\,\,(x,t)\in% \Omega_{x}\times[0,t_{f}]. 
Let f_{0} denote a given initial distribution f(x,v,0)=f_{0}(x,v) in (x,v)\in[0,1]\times\mathbb{R}. We impose periodic boundary conditions on x for the transport equation (1.1),
f(0,v,t)=f(1,v,t)\quad\forall\,\,(v,t)\in\mathbb{R}\times[0,t_{f}], 
and also for the Poisson equation (1.2); i.e.,
(1.3)  \Phi(0,t)=\Phi(1,t),\quad\forall\,\,t\in[0,t_{f}]. 
To ensure the wellposedness of the Poisson problem we add the compatibility (or normalizing) condition
(1.4)  \displaystyle{\int_{0}^{1}\rho(x,t)dx=\int_{0}^{1}\int_{\mathbb{R}}f(x,v,t)% dvdx=1},\quad\forall\,\,t\in[0,t_{f}], 
which is the condition for total charge neutrality. To guarantee the uniqueness of its solution \Phi (otherwise is determined only up to a constant), we fix the value of \Phi at a point. We set
(1.5)  \Phi(0,t)=0\quad\forall\,\,t\in[0,t_{f}]. 
Notice that (1.4) express that the total charge of the system is preserved in time.
Through the paper we are only concerned with compactly supported solutions f of problem (1.6)(1.2). We assume that a bounded set \Omega_{v}\subset\mathbb{R} such that
\mbox{supp}(f(x,v,0))\cup\mbox{supp}(f(x,v,t))\subseteq\Omega_{x}\times\Omega_% {v}\;,\qquad\forall\,t\in[0,t_{f}]\;, 
and so the Vlasov equation (1.1), can be (and will be) regarded in \Omega_{x}\times\Omega_{v}\times[0,t_{f}]:
(1.6)  \frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}\Phi_{x}\frac{% \partial f}{\partial v}=0\quad(x,v,t)\,\,\in\,\,\Omega_{x}\times\Omega_{v}% \times[0,t_{f}]\;. 
The charge density is accordingly defined by
(1.7)  \rho(x,t)=\displaystyle{\int_{\Omega_{v}}f(x,v,t)dv}\quad\forall\,\,(x,t)\in% \Omega_{x}\times[0,t_{f}]. 
We define the total energy of the system as
(1.8)  \mathcal{E}(t)=\int_{\Omega}f(x,v,t)\frac{v^{2}}{2}dxdv+\int_{\Omega_{x}}% \frac{1}{2}\Phi_{x}(x,t)^{2}dx\quad\forall\,t\in[0,t_{f}]. 
The first term in the above definition represents the kinetic energy; the second, the potential energy of the system.
1.1.1. Properties
The VlasovPoisson system preserves in time many physical observables. We now briefly revise some:

Mass conservation: as already mentioned, the total charge of the system is preserved:
(1.9) \frac{d}{dt}\int_{\Omega}f(x,v,t)\,dx\,dv=0\qquad\forall\,t\in[0,t_{f}]\;. 
L^{p}conservation: noting that {\rm div}_{x,v}\left([v,\Phi_{x}(x,t)]\right)\equiv 0 one can deduce straightaway the conservation of all L^{p}norms of the distribution function:
(1.10) \frac{d}{dt}\int_{\Omega}\f(x,v,t)\^{p}\,dx\,dv=0\qquad\forall\,t\in[0,t_{f}% ]\;. We will be particularly concerned with p=1,2.

Total Energy: Following [dolbeault], one can also show the following energy apriori estimate:
(1.11) \frac{d}{dt}\mathcal{E}(t)=\frac{d}{dt}\left(\int_{\Omega}f(x,v,t)\frac{v^{2% }}{2}\,dx\,dv+\int_{\Omega_{x}}\frac{1}{2}E(x,t)^{2}dx\right)=0\qquad\forall% \,t\in[0,t_{f}]\;, where we have already used the definition of the electrostatic field E(x,t)=\Phi_{x}(x,t) (compare with (1.8)).
1.2. Basic notation and preliminaries for the numerical methods
Let \{\mathcal{T}_{h}\} be a family of partitions of our computational/physical domain \Omega=\Omega_{x}\times\Omega_{v}=\Omega_{x}\times[L,L], which we assume to be regular [ciar2] and made of rectangles. Each cartesian mesh \mathcal{T}_{h} is defined as
\mathcal{T}_{h}:=\left\{T_{ij}=I_{i}\times J_{j},\quad 1\leq i\leq N_{x},\,\,1% \leq j\leq N_{v}\,\right\}, 
where
I_{i}=[x_{i1/2},x_{i+1/2}]\quad\forall\,i=1,\ldots,N_{x};\qquad J_{j}=[v_{j1% /2},v_{j+1/2}]\quad\forall\,j=1,\ldots,N_{v}\;. 
The mesh sizes h_{x} and h_{v} relative to the partition are defined as
0<h_{x}=\max_{1\leq i\leq N_{x}}h_{i}^{x}:=x_{i+1/2}x_{i1/2},\quad 0<h_{v}=% \max_{1\leq j\leq N_{v}}h_{j}^{v}:=v_{i+1/2}v_{i1/2}\;, 
with h_{i}^{x} and h_{j}^{v} denoting the cell lengths of I_{i}
and J_{j}, respectively. The mesh size of the partition is
defined as h=\max{(h_{x},h_{v})}. The shape regularity assumption implies that \exists\,c_{1},\,c_{2}>0 constants independent of h such that c_{1}h_{v}\leq\ h_{x}\leq c_{2}\,h_{v}.
We assume that v=0 corresponds to a node of the partition along the vaxis, i.e.,
v_{j1/2}=0 for some j in the partition of \Omega_{v}=[L,L]. We denote by \{\mathcal{I}_{h}\} the family of partitions of the interval \Omega_{x}: \mathcal{I}_{h}:=\left\{\,\,I_{i}\,:\,\,1\leq i\leq N_{x}\,\right\}.
For k\geq 1, let \mathbb{P}^{k}(I_{i}) be the space of polynomials of degree up to k, and let \mathbb{Q}^{k}(T_{ij}) be the space
of polynomials of degree at most k in each variable ((x,v)). We define the finite element spaces:
(1.12)  \displaystyle V_{h}^{k}  \displaystyle=  \displaystyle\left\{\psi\in L^{2}(\mathcal{I})\,\,:\quad\psi\in\mathbb{P}^{k}(% I_{i}),\,\,\,\forall\,I_{i}\,,\,i=1,\ldots N_{x},\right\},  
(1.13)  \displaystyle\mathcal{Z}_{h}^{k}  \displaystyle:=  \displaystyle\left\{\xi\in L^{2}(\Omega)\,\,:\quad\xi\in\mathbb{Q}^{k}(T_{ij})% ,\,\,\,\forall\,T_{ij}=I_{i}\times J_{j},\,\,\forall i\,,j\,\right\},  
(1.14)  \displaystyle W_{h}^{k}  \displaystyle=  \displaystyle\left\{\chi\in\mathcal{C}^{0}(\mathcal{I})\,\,:\quad\chi\in% \mathbb{P}^{k}(I_{i}),\,\,\,\forall\,I_{i}\,,\,i=1,\ldots N_{x},\right\}\cap L% ^{2}(\mathcal{I})/\mathbb{R}\;. 
As is usual in the DG methods, we now introduce the the trace operators. We denote by (\varphi_{h})_{i+1/2,v}^{+} and (\varphi_{h})_{i+1/2,v}^{} the values of \varphi_{h} at (x_{i+1/2},v) from the right cell I_{i+1}\times J_{j} and from the left cell I_{i}\times J_{j}, respectively;
(\varphi_{h})_{i+1/2,v}^{\pm}=\displaystyle{\lim_{\varepsilon\downarrow 0}{% \varphi_{h}(x_{i+1/2}\pm\varepsilon,v)}}\;,\quad(\varphi_{h})_{x,j+1/2}^{\pm}=% \displaystyle{\lim_{\varepsilon\downarrow 0}{\varphi_{h}(x,v_{j+1/2}\pm% \varepsilon)}}\;, 
for all (x,v)\in\mathcal{I}\times\mathcal{J} or in shorthand notation
(1.15)  (\varphi_{h})_{i+1/2,v}^{\pm}=\varphi_{h}(x^{\pm}_{i+1/2},v)\;,\qquad(\varphi_% {h})_{x,j+1/2}^{\pm}=\varphi_{h}(x,v^{\pm}_{j+1/2})\;, 
for all (x,v)\in I_{i}\times J_{j}. The jump [\![\,\cdot\,]\!] and average \{\cdot\} trace operators of \varphi_{h} at (x_{i+1/2},v),\,\,\forall\,v\in J_{j} are defined by
\displaystyle[\![\,\varphi_{h}\,]\!]_{i+1/2,v}  \displaystyle:=(\varphi_{h})_{i+1/2,v}^{+}(\varphi_{h})_{i+1/2,v}^{}  \displaystyle\forall\varphi_{h}\in\mathcal{Z}_{h}^{k}\;,  
\displaystyle\{\varphi_{h}\}_{i+1/2,v}  \displaystyle:=\frac{1}{2}\left[(\varphi_{h})_{i+1/2,v}^{+}+(\varphi_{h})_{i+1% /2,v}^{}\right]  \displaystyle\forall\varphi_{h}\in\mathcal{Z}_{h}^{k}\;. 
For k\geq 0, let P^{k}:L^{2}(\mathcal{I})\longrightarrow V_{h}^{k} be the standard L^{2} orthogonal projection onto the finite element space V_{h}^{k} defined locally, i.e., for each 1\leq i\leq N_{x},
(1.16)  \int_{I_{i}}\left(P^{k}(w)w\right)q_{h}\,dx=0\qquad\forall q_{h}\in\mathbb{P}% ^{k}(I_{i})\;. 
By definition the projection is stable in L^{2}(\mathcal{I})
(1.17)  \P^{k}(w)\_{L^{2}(\mathcal{I}_{h})}\leq\w\_{L^{2}(\mathcal{I})}\qquad% \forall\,w\in L^{2}(\mathcal{I}). 
We denote by \mathcal{P}_{h}:L^{2}(\Omega)\longrightarrow\mathcal{Z}_{h}^{k} the corresponding two dimensional L^{2}orthogonal projection; defined by \mathcal{P}_{h}(w)=(P^{k}_{x}\otimes P^{k}_{v})(w); i.e., for all i and j,
(1.18)  \int_{I_{i}}\int_{J_{j}}\left(\mathcal{P}_{h}(w(x,v))w(x,v)\right)\varphi_{h}% (x,v)\,dv\,dx=0\quad\forall\varphi_{h}\in\mathbb{P}^{k}(I_{i})\otimes\mathbb{P% }^{k}(J_{j})\;. 
Also from its definition, its L^{2}stability follows immediately.
2. Discontinuous Galerkin methods for the VlasovPoisson system: semidiscrete methods
In this section, we introduce the DG methods we consider for approximating the VlasovPoisson system. The methods are those proposed in [acs0, acs1], but with some slight variation required for practical computations. Following [acs0, acs1] we first describe the schemes for the Vlasov equation, proposing several options to modify the methods in [acs0] so that they allow for a feasible implementation. We then discuss the approximation of the Poisson problem (again following closely [acs0]). We close the section by discussing the main properties of the introduced schemes. Throughout the whole section we focus on the space discretization.
2.1. Discontinuous Galerkin approximation to the Vlasov equation
We now describe the DG methods we consider to approximate the Vlasov equation (1.6).
For the time being, we assume we are given a finite element (conforming or nonconforming) approximation of degree r to the electrostatic field E(x,t)=\Phi_{x}(x,t), which we denote by E_{h}\in\mathcal{W}_{h}. By E_{h}^{i} we refer to its restriction to I_{i}. The properties and characterization of E_{h} are discussed in next subsection.
We denote by f_{h}(0)=\mathcal{P}_{h}(f(x,v,0)) the approximation to the initial data f(x,v,0) computed using the orthogonal L^{2}projection onto the space \mathcal{Z}_{h}^{k}. Since the Vlasov equation is a transport equation, we construct the DG method in the usual way: given E_{h}\in\mathcal{W}_{h} find f_{h}:[0,t_{f}]\longrightarrow\mathcal{Z}_{h}^{k} such that
(2.1)  \displaystyle\sum_{i=1}^{N_{x}}\displaystyle\sum_{j=1}^{N_{v}}\mathcal{B}^{h}_% {ij}(E_{h};f_{h},\varphi_{h})=0\quad\forall\varphi_{h}\in\mathcal{Z}_{h}^{k}\;, 
where the bilinear form \mathcal{B}^{h}_{ij}(E_{h};f_{h},\varphi_{h}) is defined for each i, j and \varphi_{h}\in\mathcal{Z}_{h}^{k} as:
\displaystyle\mathcal{B}_{ij}(E_{h};f_{h},\varphi_{h})=  \displaystyle\displaystyle{\int_{T_{ij}}\frac{\partial f_{h}}{\partial t}% \varphi_{h}\,dv\,dx\int_{T_{ij}}vf_{h}\frac{\partial\varphi_{h}}{\partial x}% \,dv\,dx+\int_{T_{ij}}E_{h}^{i}f_{h}\frac{\partial\varphi_{h}}{\partial v}\,dv% \,dx}  
(2.2)  \displaystyle\displaystyle{+\int_{J_{j}}\left[(\widehat{(vf_{h})}\varphi^{}_{% h})_{i+1/2,v}(\widehat{(vf_{h})}\varphi^{+}_{h})_{i1/2,v}\right]dv}  
\displaystyle\displaystyle{\int_{I_{i}}\left[\left(\widehat{\left(E^{i}_{h}f_% {h}\right)}\varphi^{}_{h}\right)_{x,j+1/2}\left(\widehat{\left(E^{i}_{h}f_{h% }\right)}\varphi^{+}_{h}\right)_{x,j1/2}\right].dx}, 
In (2.2) we have used the short hand notation given in (1.15). The numerical fluxes are defined using the upwind flux:
(2.3)  \displaystyle\widehat{vf_{h}}  \displaystyle=\left\{\begin{array}[]{cc}v\,f_{h}^{}&\mbox{ if }v\geq 0\\ v\,f_{h}^{+}&\mbox{ if }v<0\end{array}\right.  \displaystyle\widehat{vf_{h}}=\{vf_{h}\}\frac{v}{2}[\![\,f_{h}\,]\!]\;,  
(2.4)  \displaystyle\widehat{E_{h}^{i}f_{h}}  \displaystyle=\left\{\begin{array}[]{cc}E_{h}^{i}\,f_{h}^{+}&\mbox{ if }% \mathcal{P}^{0}(E_{h}^{i})\geq 0\\ E_{h}^{i}\,f_{h}^{}&\mbox{ if }\mathcal{P}^{0}(E_{h}^{i})<0\end{array}\right.  \displaystyle\widehat{E_{h}^{i}f_{h}}=\{E_{h}^{i}f_{h}\}+\mbox{sign}\left(% \mathcal{P}^{0}(E_{h}^{i})\right)\cdot\frac{E_{h}^{i}}{2}[\![\,f_{h}\,]\!]\;. 
At the boundary \partial\Omega, the numerical fluxes are taken as
(\widehat{vf_{h}})_{1/2,v}=(\widehat{vf_{h}})_{N_{x}+1/2,v},\quad(\widehat{E_{% h}^{i}f_{h}})_{x,1/2}=(\widehat{E_{h}^{i}f_{h}})_{x,N_{v}+1/2}=0,\,\forall\,(x% ,v)\in\mathcal{I}\times\mathcal{J}, 
so that the periodicity in x and the compactness in v are
reflected.
Note that the numerical fluxes as defined in (2.3) and (2.4) are consistent.
Observe that, unlike in [acs0, acs1] the definition (2.4) of the upwind flux \widehat{E_{h}^{i}f_{h}} involves a condition on the sign(\mathcal{P}^{0}(E_{h}^{i})), rather than on sign(E_{h}^{i}). Obviously if E_{h}^{i} does not vanish inside I_{i} (and so it does not change sign inside I_{i}), sign(E_{h}^{i})=\mbox{sign}\left(\mathcal{P}^{0}(E_{h}^{i})\right) and the classical definition of the upwind flux is recovered. However, since E_{h} is a piecewise polynomial of degree k+1 approximation to the electrostatic field, it will in general change sign in some elements I_{i} of the partition \mathcal{I}_{h}. The classical definition would require to construct (at each time step) a partition of \Omega_{x} that adapts to the changes of sign of E_{h} by locating the zeros of E_{h} at nodes of the desired partition. Such process, although feasible in one dimension, might become too expensive and complicate unnecessarily the whole solution method for the ValsovPoisson system (specially in higher dimensions).
The definition (2.4) is considered for computational purposes. It allows to avoid
computing the zeros of E_{h} and remeshing, at each time step, the partition \mathcal{I}_{h}.
For our computations of the onedimensional problem we have also examined two other variants of the numerical flux \widehat{E_{h}^{i}f_{h}} defined in (2.4), that do not require remeshing, although they require a control on the sign of E_{h}. We also show how this control on the sign(E_{h}) can be done efficiently. The first variant we consider is given by
(2.5)  \widehat{E_{h}^{i}f_{h}}=\left\{\begin{aligned} \displaystyle E_{h}^{i}\,f_{h}% ^{+}&\displaystyle\qquad\mbox{ if }E_{h}^{i}>0&&\\ \displaystyle E_{h}^{i}\,f_{h}^{}&\displaystyle\qquad\mbox{ if }E_{h}^{i}<0&&% \\ \displaystyle\{E_{h}^{i}\,f_{h}\}_{\omega}&\displaystyle\qquad\mbox{ if }% \exists\,x^{\ast}\in I_{i}\quad\mbox{such that}\quad E_{h}^{i}(x^{\ast})=0&&% \end{aligned}\right. 
where we have used the weighed average
(2.6)  \{E_{h}^{i}\,f_{h}\}_{\omega}=\omega^{+}E_{h}^{i}f_{h}^{+}+\omega^{}E_{h}^{i}% f_{h}^{}\qquad\omega^{+}\;,\,\omega^{}\in[0,1]\quad\omega^{+}+\omega^{}=1. 
The parameter \omega should be chosen so that the amount of upwind is tuned. Although based on heuristics, in our computations we have found that a good choice is given by
(2.7)  \omega^{+}=\frac{\max_{I_{i}}E_{h}}{\max_{I_{i}}E_{h}+\min_{I_{i}}E_{h}}% \;,\qquad\omega^{}=\frac{\min_{I_{i}}E_{h}}{\max_{I_{i}}E_{h}+\min_{I_{i% }}E_{h}}\;. 
The definition of the numerical flux (2.5) can be rewritten in the compact form:
(2.8)  \widehat{E_{h}^{i}f_{h}}=\left\{\begin{aligned} \displaystyle\{E_{h}^{i}f_{h}% \}+\frac{E_{h}^{i}}{2}[\![\,f_{h}\,]\!]&\displaystyle\qquad\mbox{ if }% \nexists\,x^{\ast}\in I_{i}\quad\mbox{such that}\quad E_{h}^{i}(x^{\ast})=0&&% \\ \displaystyle\{E_{h}^{i}f_{h}\}+E_{h}^{i}(\omega^{+}\frac{1}{2})[\![\,f_{h}\,% ]\!]&\displaystyle\qquad\mbox{ if }\exists\,x^{\ast}\in I_{i}\quad\mbox{such % that}\quad E_{h}^{i}(x^{\ast})=0&&\end{aligned}\right. 
Observe that this definition of the numerical flux is also consistent.
The last variant we consider is defined by:
(2.9)  \widehat{E_{h}^{i}f_{h}}=\left\{\begin{array}[]{cc}E_{h}^{i}\,f_{h}^{+}&\mbox{% if }E_{h}^{i}>0\\ E_{h}^{i}\,f_{h}^{}&\mbox{ if }E_{h}^{i}<0\\ \mathcal{P}^{0}(E_{h}^{i})f_{h}^{+}&\mbox{ if }\mathcal{P}^{0}(E_{h}^{i})>0% \mbox{ and }\exists\,x^{\ast}\in I_{i}\quad\mbox{such that}\quad E_{h}^{i}(x^{% \ast})=0\\ \mathcal{P}^{0}(E_{h}^{i})f_{h}^{}&\mbox{ if }\mathcal{P}^{0}(E_{h}^{i})<0% \mbox{ and }\exists\,x^{\ast}\in I_{i}\quad\mbox{such that}\quad E_{h}^{i}(x^{% \ast})=0.\\ \end{array}\right. 
Note that the numerical flux defined above in (2.9) is not consistent (it fails to be consistent in those elements where E_{h}^{i} changes sign, where we commit an error of order O(h)).
Notice that both the weighted average approach (2.5) and the last definition (2.9) require the knowledge of those elements of the partition \mathcal{I}_{h} where E_{h} vanishes. This information can be obtained very easily (at least in one dimension), by checking the sign of coefficients of E_{h} expanded in a basis with Bernstein polynomials^{1}^{1}1Bernstein polynomials are nonnegative at everypoint of their domain. See the Appendix LABEL:ap0 for further details on Bernstein polynomials and how the detection of change of sign is implemented.
Mimicking (1.7), we define the discrete density, \rho_{h}(x,t):
(2.10)  \rho_{h}(x,t)=\int_{\mathcal{J}}f_{h}(x,v,t)\,dv=\displaystyle{\sum_{j}\int_{J% _{j}}f_{h}(x,v,t)\,dv}\quad\forall\,\,x\,\in\,\mathcal{I},\quad\forall\,\,t\,% \in\,\,[0,t_{f}]. 
One of the nice properties of DG schemes, is that the conservation of the total mass is satisfied by construction. Next Lemma guarantees that the DG scheme (2.1)(2.2) with fluxes (2.3) and either (2.4) or (2.5) preserves the total charge:
Lemma 2.1.
Particle or Mass Conservation: Let f_{h}\in\mathcal{C}^{1}([0,t_{f}];\mathcal{Z}_{h}^{k}), with k\geq 0, be the DG approximation to f, satisfying (2.1)(2.2), with numerical fluxes defined as in (2.3) and either (2.4) or (2.5) or (2.9). Then, for all t\in[0,t_{f}],
(2.11)  \displaystyle\sum_{i,j}\int_{T_{ij}}f_{h}(t)\,dv\,dx=\displaystyle\sum_{i,j}% \int_{T_{ij}}f_{h}(0)\,dv\,dx=\displaystyle\sum_{i,j}\int_{T_{ij}}f_{0}\,dv\,% dx=1. 
Although standard, we provide here the proof of the above Lemma for the sake of completeness.
Proof.
(Proof of Lemma 2.1). Since f_{h}(0)=\mathcal{P}_{h}(f_{0}) it follows from the mass conservation of the continuous VP system (1.4) and the definition of the L^{2}projection (1.18) that
(2.12)  \displaystyle\sum_{i,j}\int_{T_{ij}}f_{h}(0)\,dv\,dx=\displaystyle\sum_{i,j}% \int_{T_{ij}}\mathcal{P}_{h}(f_{0})\,dv\,dx=\displaystyle\sum_{i,j}\int_{T_{ij% }}f_{0}\,dv\,dx=1. 
Now, let T_{ij} be any arbitrary but fixed element in \mathcal{T}_{h}. By setting in (2.2) \varphi_{h}=1 in T_{ij} and \varphi_{h}=0 elsewhere we find,
\displaystyle\mathcal{B}_{ij}(E_{h}{\color[rgb]{1,0,0};}f_{h},1)=  \displaystyle\,\frac{d}{dt}\int_{T_{ij}}f_{h}\,dv\,dx+\int_{J_{j}}[\widehat{(% vf_{h})}_{i+1/2,v}\widehat{(vf_{h})}_{i1/2,v}]\,dv  
\displaystyle\int_{I_{i}}[\widehat{(E^{i}_{h}f_{h})}_{x,j+1/2}\widehat{(E^{i% }_{h}f_{h})}_{x,j1/2}]\,dx\;, 
where we have already used that such \varphi_{h} obviously satisfies (\varphi_{h})_{i+1/2,v}^{}=(\varphi_{h})_{i1/2,v}^{+}=1 at the boundaries of T_{ij}. Now, the above equation obviously holds for any i,j, since the choice of T_{ij} was arbitrary. Therefore summing over all i and j the above equation, the flux terms telescope and there is no boundary term left because of the periodic (for i) and compactly supported (for j) boundary conditions. Substitution now in (2.1) gives
0=\displaystyle\sum_{i,j}\mathcal{B}_{ij}(E_{h};f_{h},1)=\frac{d}{dt}% \displaystyle\sum_{i,j}\int_{T_{ij}}f_{h}\,dv\,dx=0, 
and so integrating in time and using (2.12) we reach (2.11). ∎
2.2. Finite element approximation of the electrostatic field E
We now describe the methods we consider for approximating the electrostatic field E(x,t)=\Phi_{x}(x,t). The discrete Poisson problem reads:
(2.13)  (\Phi_{h})_{xx}=1\rho_{h}\quad x\in\Omega_{x},\qquad\Phi_{h}(1,t)=\Phi_{h}(0,% t). 
The well posedness of the above discrete problem is guaranteed by (2.11) from Lemma 2.1 which in particular implies
(2.14)  (\Phi_{h})_{x}(1,t)=(\Phi_{h})_{x}(0,t). 
To ensure the uniqueness of the solution we also set \Phi_{h}(0,t)=0.
Since in the Vlasov equation, the transport depends on E, to approximate (2.13) we consider a mixed finite element approach. For that purpose, we first rewrite problem (2.13) as a first order system:
(2.15)  E_{h}=\frac{\partial\Phi_{h}}{\partial x}\quad x\in\Omega_{x};\qquad\frac{% \partial E_{h}}{\partial x}=\rho_{h}1\quad x\in\Omega_{x}\, 
with boundary condition \Phi_{h}(0,t)=\Phi_{h}(1,t)=0. We consider the following methods:
2.2.1. Mixed Finite element approximation:
we consider the onedimensional version of RaviartThomas elements, RT{}_{k}\,\,k\geq 1 [ravtom0, brezzifortin]. In 1D the mixed finite element spaces turn out to be the (W_{h}^{k+1},V_{h}^{k})finite element spaces. Note that in particular, \frac{d}{dx}(W_{h}^{k+1})=V_{h}^{k}. For k\geq 0 the scheme reads: find (E_{h},\Phi_{h})\in W_{h}^{k+1}\times V_{h}^{k} such that
(2.16)  \displaystyle\int_{\mathcal{I}}E_{h}\,z\,dx+\int_{\mathcal{I}}\Phi_{h}\,z_{x}% \,dx=0  \displaystyle\forall\,z\in W_{h}^{k+1},  
(2.17)  \displaystyle\int_{\mathcal{I}}(E_{h})_{x}\,p\,dx=\int_{\mathcal{I}}(\rho_{h}% 1)p\,dx  \displaystyle\quad\forall\,p\in V_{h}^{k}. 
We also refer to [babuskanarashi], where the lowest order case was studied for the onedimensional Poisson problem.
2.2.2. Local Discontinuous Galerkin (LDG) method:
The DG approximation to the first order system (2.15) reads: find (E_{h},\Phi_{h})\in V_{h}^{k+1}\times V_{h}^{k+1} such that for all i:
(2.18)  \displaystyle\int_{I_{i}}E_{h}z\,dx=\int_{I_{i}}\Phi_{h}z_{x}\,dx+[(\widehat{% \Phi_{h}}z^{})_{i+1/2}(\widehat{\Phi_{h}}z^{+})_{i1/2}]  \displaystyle\forall\,z\in V_{h}^{k+1},  
(2.19)  \displaystyle\int_{I_{i}}E_{h}p_{x}\,dx\left[(\widehat{E_{h}}p^{})_{i+1/2}(% \widehat{E_{h}}p^{+})_{i1/2}\right]=\int_{I_{i}}(\rho_{h}1)p\,dx  \displaystyle\forall\,p\in V_{h}^{k+1}\;. 
The numerical fluxes (\widehat{\Phi_{h}})_{i1/2} and (\widehat{E_{h}})_{i1/2} for the LDG method are defined by:
(2.20)  \left\{\begin{aligned} \displaystyle(\widehat{\Phi_{h}})_{i1/2}&\displaystyle% =\{\Phi_{h}\}_{i1/2}c_{12}[\![\,\Phi_{h}\,]\!]_{i1/2}\;,&&\\ \displaystyle(\widehat{E_{h}})_{i1/2}&\displaystyle=\{E_{h}\}_{i1/2}+c_{12}[% \![\,E_{h}\,]\!]_{i1/2}+c_{11}[\![\,\Phi_{h}\,]\!]_{i1/2}\;,&&\end{aligned}\right. 
where c_{11}=c\,(k+1)^{2}h_{x}^{1} and c_{12}=1/2. At the boundary nodes due to periodicity in x we impose
(\widehat{\Phi_{h}})_{1/2}=(\widehat{\Phi_{h}})_{N_{x}+1/2},\quad(\widehat{E_{% h}})_{1/2}=(\widehat{E_{h}})_{N_{x}+1/2}. 
The method was first introduced in [cockshu98] for a time dependent convection diffusion problem with c_{11}=O(1). For the Poisson problem it has been considered in [faithcock] in the onedimensional case, and in [ilaria0] for higher dimensions.