A WENObased Method of Line Transpose Approach for Vlasov Simulations
Andrew Christlieb^{1}^{1}1 Department of Mathematics and Electrical Engineering, Michigan State University, East Lansing, MI, 48824. Email: andrewch@math.msu.edu , Wei Guo ^{2}^{2}2 Department of Mathematics, Michigan State University, East Lansing, MI, 48824. Email: wguo@math.msu.edu, Yan Jiang^{3}^{3}3Department of Mathematics, Michigan State University, East Lansing, MI, 48824. Email: jiangyan@math.msu.edu
Abstract.
In this paper, a high order implicit Method of Line Transpose (MOL) method based on a weighted essentially nonoscillatory (WENO) methodology is developed for onedimensional linear transport equations and further applied to the VlasovPoisson (VP) simulations via dimensional splitting. In the MOL framework, the time variable is first discretized by a diagonally implicit strongstabilitypreserving RungeKutta method, resulting in a boundary value problem (BVP) at the discrete time levels. Then an integral formulation coupled with a high order WENO methodology is employed to solve the BVP. As a result, the proposed scheme is high order accurate in both space and time and free of oscillations even though the solution is discontinuous or has sharp gradients. Moreover, the scheme is able to take larger time step evolution compared with an explicit MOL WENO scheme with the same order of accuracy. The desired positivitypreserving (PP) property of the scheme is further attained by incorporating a newly proposed high order PP limiter. We perform numerical experiments on several benchmarks including linear advection, solid body rotation problem; and on the Landau damping, twostream instabilities, bumpontail, and plasma sheath by solving the VP system. The efficacy and efficiency of the proposed scheme is numerically verified.
Key Words: Method of Line Transpose, implicit RungeKutta methods, weighted essentially nonoscillatory methodology, high order accuracy, VlasovPoisson, demensional splitting, positivitypreserving
1 Introduction
Transport modeling plays an important role in a wide range of applications, including fluid dynamics, plasma physics, atmospheric sciences, among many others. The main focus of this paper is to develop of efficient high order transport solvers for the Vlasov simulations. As a highdimensional kinetic transport equation, the VlasovPoisson (VP) system is considered a fundamental model in plasma physics which describes the dynamics of charged particles due to the selfconsistent electric force:
(1.1a)  
(1.1b) 
where describes the probability of finding a particle with velocity at position at time , is the electrostatic field, is the selfconsistent electrostatic potential, and is the electron charge density and the 1 represents the uniformly distributed infinitely massive ions in the background.
The VP simulation is challenging. Besides the famous curse of dimensionality, the underlying system may develop filamentation solution structures, since the longrange forces drive the solution out of thermoequilibrium. Hence, it is of importance that the used numerical scheme is able to effectively capture filamentation structures without producing spurious oscillations. Further, the VP system conserves several physical quantities including the total mass, momentum, and energy, which requires that the numerical scheme converses those invariants on the discrete level. In addition, as a probability density function in phase space, the unknown function always remains nonnegative. In practice, the positivity preserving property is another highly desired property for the design of Vlasov solver, since negative solutions may trigger unphysical phenomenon in plasmas.
There exist a large amount of successful VP solvers in the literature, such as the prominent particleincell (PIC) method [4, 23, 41] , which evolves the solution by following the trajectories of some sampled macroparticles, Eulerian type approaches [18, 50, 12, 11], for which the state variable is advanced according to the underlying equations on a fixed numerical grid, and semiLagrangian approach [9, 29, 2, 39, 19, 3, 17, 16, 34, 30, 15, 22], which updates the solution by tracking the characteristics on a fixed mesh grid. However, these popular schemes still have limitations. For example, the PIC method is known to suffer from low order sampling noise; Eulerian approaches are subject to stringent time step restriction associated with explicit time stepping; general boundary conditions are difficult to handle in a semiLagrangian setting. Also note that many of these existing methods are in the Method of Line (MOL) framework, meaning that the spatial variable is first discretized, then the numerical solution is updated in time by coupling a suitable time integrator. In this paper, we consider an alternative approach to advance the solution, which is called the Method of Line Transpose (MOL) or Rothes’s method in the literature [35, 36]. As its name implies, such an approach is defined in an orthogonal fashion of the MOL schemes, i.e., the discretization is first carried out for temporal variable, resulting in a boundary value problem (BVP) at the discrete time levels. Then a preferred BVP solver can be further applied to advance the solution. For example, in [6], an implicit integral formulation based on the MOL framework for onedimensional wave equations is developed and analyzed. It is proven that the scheme is Astable, hence allowing for extra large time step evolution. In [8], the computational cost of the scheme is reduced from to , where is the number of discrete mesh points, by taking advantage of an important property of the analytical solution of the wave equations. More recently, a successive convolution technique [5] is developed to enhance the temporal accuracy to arbitrary high order. A notable advantage of such a scheme is that, even though it is implicit in time, we do not need to explicitly solve a linear system, while the BVP is inverted analytically in an integral formulation. The method is further extended to multidimensional wave equations by utilizing the alternating direction implicit strategy at the expense of introducing splitting errors [6]. However, the simplicity of the scheme is preserved. The recent applications of this MOL approach include the heat equation, the AllenCahn equation [7], and Maxwell’s equations [10].
In this paper, we consider formulating a numerical scheme suited for the Vlasov equation in the MOL framework in order to take advantage of its good properties such as the simplicity and the ability to take large time steps for evolution. In order to accommodate the MOT approach, the nonlinear multidimensional Vlasov equation is first split into several lowerdimensional linear transport equations. Such an idea has been widely used in semiLagrangian Vlasov simulations [9]. As we mentioned earlier, the extra difficulty of VP simulations lies in the distinctive filamentation solution structures, which requires that the employed numerical scheme is able to effectively capture sharp gradient of the solution without producing spurious oscillations even on an underresolved mesh. Note that we cannot directly apply the existing MOL schemes to the VP system, since they are all linear schemes and hence exhibit oscillations when the solution is discontinuous or has a sharp gradient. In order to control undesired oscillations, we propose to further couple two numerical ingredients. Firstly, in the MOL framework, the time variable is discretized via an implicit high order strongstabilitypreserving (SSP) RungeKutta (RK) method [21]. Secondly, a robust weighted essential nonoscillatory (WENO) methodology is incorporated when inverting the BVP. We expect that the resulting scheme is efficient, high order accurate, robust, and able to take large time step evolution (larger than the explicit counterpart). Furthermore, a high order positivitypreserving (PP) limiter suited for the proposed MOL scheme is developed to ensure such a high desired property for the numerical solution.
We also would like to point out that the proposed MOL WENO scheme shares some similarity with the semiLagrangian WENO scheme developed in [31, 32], in the sense that both schemes are high order accurate and able to effectively capture the filamentation structures with large time step. On the other hand, compared with semiLagrangian WENO schemes, the proposed method is able to conveniently handle general boundary conditions, which is not trivial for semiLagrangian approaches.
The rest of the paper is organized as follows. In Section 2, we formulate the numerical scheme for onedimensional linear transport equations by coupling several main numerical ingredients such as the MOL approach, the WENO methodology and the high order PP limiter. In Section 3, the method is extended to twodimensional transport equations via dimensional splitting. We present several benchmark tests in Section 4 to verify the performance of the proposed scheme, including linear transport, rigid body rotation, as well as Landau damping, twostream instabilities, bumpontail and plasma sheath for VP simulations. We conclude the paper in Section 5.
2 Formulation of the Scheme
We start with the following onedimensional advection equation
(2.2) 
where is constant and denotes the wave propagation speed. The boundary condition is assumed to be periodic, Dirichlet boundary condition
or Neumann boundary condition
2.1 The method of line transpose approach
As discussed in the introduction, in a MOL framework, we first discretize the time variable while leaving the spatial variable continuous. For example, we can apply the first order backward Euler method to (2.2), and obtain
i.e.,
(2.3) 
where denotes the time step and . Note that (2.3) is an initial value problem (IVP), for which we can obtain an explicit representation for solution at time level . If we assume , can be represented by
(2.4) 
where
(2.5) 
and is a constant and should be determined by the boundary condition. The superscript indicates that the characteristics traverse from the left to the right. If , then the characteristics traverse from the right to the left. Due to the symmetry, we can advance the solution as
(2.6) 
when the first order backward Euler scheme is used. Here,
(2.7) 
and again the constant is determined by the boundary condition.
We need to further discretize the spatial variable based on the semidiscrete schemes (2.4) and (2.6). Following the idea in [6, 8], the domain is discretized by uniformly distributed grid points
and denotes the mesh size. The solution is accordingly discretized at points , and we use to denote the numerical solution at location at time level . The main part of the algorithm is to evaluate the convolution integral in (2.4), (2.6) or in (2.5), (2.7) based on the point values .
We first consider the case of . Note that,
(2.8) 
where
(2.9) 
and is a positive constant. Therefore, in order to obtain denoted by , we only need to approximate the increment denoted by , and obtain via a recursive relation, i.e.,
(2.10) 
For the case of , again thanks to the symmetry, we have the following identity
(2.11) 
where
(2.12) 
and hence a similar recursive relation can be obtained
(2.13) 
where and denote the approximate values of and , respectively.
After computing the contributions from the left or right characteristic, we further need to find the global coefficient in (2.4) or in (2.6), which is used to enforce the boundary condition. Since , for a Dirichlet boundary condition, the global coefficient can be obtained by
(2.14) 
Similarly, for a Neumann boundary condition, because and , we can get the coefficients by
(2.15) 
For the periodic boundary condition, the coefficient can be determined by using the mass conservation property of the solution, i.e,
On the discrete level, we enforce
or
Then the coefficient is obtained
(2.16) 
or
(2.17) 
2.2 WENO methodology
The existing MOL schemes mainly use linear interpolation methods to compute increment and such as [6], which work well for smooth problems. However, the schemes may generate spurious oscillations when the solution has discontinuities because of the famous Gibbs phenomenon. In order to control such undesired oscillations, we will incorporate a WENO methodology in the MOL framework to evaluate and below.
The first WENO scheme was introduced by Liu et al. [27], as a third order accurate finite volume scheme to solve hyperbolic conservation laws in one dimension. Later, in [24], a general framework for arbitrary order accurate finite difference WENO schemes was designed, which is more efficient for multidimensional calculations. The schems have the advantage of attaining high order accuracy in smooth regions of the solution while maintaining sharp and essentially nonoscillatory transitions of discontinuities. More information about WENO can be found in [38].
In our work, a WENO integration procedure is used, based on the idea proposed in [13, 14, 28]. Below, a th order approximation of is provided as an example, where is an integer. The corresponding coefficients of the third order and fifth order schemes are given in the Appendix. The whole procedure consists of four steps.

On each small stencil , , which contains and , there is a unique polynomial of degree at most which interpolates at the nodes in . Then we can get the integral
(2.18) where the coefficients depend on and the cell size , but not on .

Similarly, on the big stencil , there is a unique polynomial of degree at most interpolating at the nodes in . Then we have
(2.19) and linear weights , such that
(2.20) where .

Change the linear weights into nonlinear weights , using
(2.21) with
Here, is a small number to avoid the denominator becoming zero, and we take in our numerical tests. The smoothness indicator , measuring the relative smoothness of the function in the stencil , is defined by
(2.22) 
Lastly, we have
(2.23)
The process to obtain is mirrorsymmetric to that for , with respect to the point .
When dealing with a Dirichlet or a Neumann boundary condition, we further need to extrapolate point values at the ghost points and using the point values in the domain, where
Following the idea in [40], a WENO extrapolation strategy is designed to get more robust results. For instance, when approximating the left ghost point , we consider the big stencil :

On each small stencil , , we can construct a polynomial of degree at most interpolating at the nodes in . And can be extrapolated by
(2.24) where , , , , and .

Change the weights into nonlinear weights
(2.25) with
Again, we take in our numerical tests to avoid the denominator becoming zero. And the smoothness indicator is defined by

Finally, we have
Similar, we can get the WENO extrapolation on right ghost points based on .
2.3 High order time discretization
Note that scheme (2.3) is first order accurate in time. Below, we will extend our schemes to high order accuracy via the implicit SSP RK method [21, 25].
To approximate the solutions of ordinary differential equation
(2.26) 
an stage and th order RK method, denoted as , is defined as
(2.27a)  
(2.27b) 
with the standard assumption
Denote the matrix and the vector . In our scheme, we will use the diagonally implicit (DIRK) methods, for which A is lower triangular. Denote , , , and . Then, we can get . Hence
(2.28a)  
(2.28b) 
Note that A is lower triangular. So is also lower triangular, with all diagnal elements equal to zero.
For problem (2.2), . Corresponding, each stage of (2.28a) can be rewritten as
(2.29) 
where . We can solve (2.29) with MOL method mentioned above. And finally, is obtained by (2.28b).
When problem (2.2) combining with the periodic boundary condition, the global coefficients or on each stage can be obtained by (2.16) or (2.17). For a Dirichlet boundary condition, the scheme may lose accuracy if we define or . To avoid order reduction, we follow the idea in [1]. In particular, denote for . To achieve th order of accuracy, boundary values at intermediate stages are defined as
(2.30) 
for , where
or
(2.31) 
for , where
The similar idea is used when dealing with a Neumann boundary condition.
In our numerical tests, we will use and schemes to obtain third order accuracy and fourth order accuracy in time, respectively. Corresponding coefficients are listed as appendix.
2.4 Positivitypreserving limiter
The exact solution of (2.2) satisfies the positivity preserving principle, meaning that if the initial condition , then the exact solution for any time . In practice, it is highly desirable to have numerical schemes also satisfy such a property on the discrete level. It is even more crucial for the Vlasov simulation, since negative values may trigger some nonphysical precess in plasmas [33]. Recently, Zhang el at. developed a high order positivitypreserving limiters for finite volume WENO methods and finite element DG methods, which is able to maintain high order accuracy of the original scheme and easy to implement [46, 47, 48, 49]. Later, Xu et al. [45, 26, 42, 43, 44] developed a parameter maximum principle preserving flux limiter for finite difference WENO schemes by limiting the high order flux towards a first order monotone flux. However, both limiters are exclusively designed for MOL type methods and hence cannot be directly applied in the MOL setting. Below, we will formulate a new positivitypreserving (PP) limiter suited for the proposed WENO MOL scheme, and also show that such a limiter is able to preserve positive while maintaining the original high order accuracy.
Assume that we have for all at time . After obtaining the numerical solutions , we can rewrite our scheme as a conservative form
in which the artificial numerical fluxes are recursively defined by
If we are able to modify the numerical flux into , such that
(2.32) 
then the modified solution at time level will be nonnegative automatically. Roughly speaking, the main idea of the limiter is to enforce the negative solution value back to zero by decreasing the outflow flux. For instance, when , we will reset the numerical flux if . Details of the limiter are shown below.
Firstly, we consider the case and a Dirichlet boundary condition. Numerical flux is modified into , starting from to ,
(2.33) 
where
and , due to the Dirichlet boundary condition. For a Nuemann boundary condition, we can also use (2.33) to modify fluxes. Once the modified fluxes are obtained, the positivitypreserving solution is computed by (2.32).
For periodic boundary conditions, however, we need an additional step in order to preserve the positivity and periodicity of the solution. In particular, we first modify the numerical flux into by (2.33) i.e., starting from into , with
where
For periodicity, we further set . Note that since the flux at is modified, the solution at may become negative. Hence, we need screen the solution one more time by the following procedure. Starting , we compute
If we let
Otherwise, let and for any with , let Finally, we are able to obtain the positivitypreserving solution via (2.32). Note that must be less than and hence .
When , wind blows from right to left, and we will reset the numerical flux if . Similar idea can be used for different boundary conditions.
Below, we show that the PPlimiter can preserve positivity while maintaining the previous high order. In particular, we have the following proposition.
Proposition 2.1.
Proof. Without loss of generality, we only consider (2.2) with and a Dirichlet boundary condition is imposed for brevity. Let be the first solution value that is less than zero starting from , i.e., and Then we have and such that Recalling that and and hence Also note that the modification of flux
is a negative high order quantity.
Now, we consider the effect of the modification of flux to the next point value in terms of accuracy. First, we need to compute
Several cases need to be taken into account.
Case (2) If Then, by (2.33), and . There are two cases we need to treat differently. If , then by the same argument we used to treat the case at , we have and hence . If , then we have
Further, since and , combining the above estimate, we have
Consequently,
We also have the following modification of the flux
(2.34) 
For both cases, we have shown that i.e., the PPlimiter preserves the order of accuracy. Moreover, for case (1) since the flux is not changed, we can move on to consider the solution value at with the same procedure as the case at . The situation of case (2) is more complicated when considering the subsequent effect of the modification of flux (2.34). Again, we need to consider two cases.
Case (2.1) If Similar to case (1) above, the flux at should remain the same, i.e., . By (2.34)
and consequently,
Case (2.2) If . We can repeat the analysis of case (2) and obtain that
which gives the estimate
and the modification of the flux
(2.35) 
More generally, assume that the fluxes are successively modified with at most constant many times, i.e., there exists a constant independent of and , such that if the fluxes for are successively changed based on