Positivity preserving semi-Lagrangian ... - Semantic Scholar

Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system

Jing-Mei Qiu1 and Chi-Wang Shu

2

Abstract Semi-Lagrangian (SL) methods have been very popular in the Vlasov simulation community [28, 2, 3, 18, 31, 4, 23, 25]. In this paper, we propose a new Strang split SL discontinuous Galerkin (DG) method for solving the Vlasov equation. Specifically, we apply the Strang splitting for the Vlasov equation [5], as a way to decouple the nonlinear Vlasov system into a sequence of linear equations. To evolve the decoupled linear equations, we propose to couple the SL framework with the semi-discrete DG formulation. The proposed SL DG method is free of time step restriction compared with the Runge-Kutta DG method, which is known to suffer from numerical time step limitation with relatively small CFL numbers according to linear stability analysis. We apply the recently developed positivity preserving (PP) limiter [36], which is a low-cost black box procedure, to our scheme to ensure the positivity of the unknown probability density function without affecting the high order accuracy of the base SL DG scheme. We analyze the stability and accuracy properties of the SL DG scheme by establishing its connection with the direct and weak formulations of the characteristics/Lagrangian Galerkin method [22]. The quality of the proposed method is demonstrated via basic test problems, such as linear advection and rigid body rotation, and via classical plasma problems, such as Landau damping and the two stream instability.

Keywords: Vlasov simulation, discontinuous Galerkin method, semi-Lagrangian method, positivity preserving method, characteristic/Lagrangian Galerkin method. 1

Department of Mathematical and Computer Science, Colorado School of Mines, Golden, 80401. E-mail: [email protected]. Research supported by NSF grant DMS-0914852. 2 Division of Applied Mathematics, Brown University, Providence, 02912. E-mail: [email protected]. Research supported by AFOSR grant FA9550-09-1-0126 and NSF grant DMS-0809086.

1

1

Introduction

In this paper, we propose a semi-Lagrangian (SL) discontinuous Galerkin (DG) method for the 1-D advection equation with general variable coefficient ut + (a(x, t)u)x = 0,

(1.1)

and further extend it to the high dimensional problem ut + ∇x (a(x, t)u) = 0,

(1.2)

via Strang splitting. The design of the proposed method is motivated by solving the multidimensional Vlasov-Poisson (VP) system via Strang split SL approaches originally proposed in [5], and further developed in the finite volume and finite difference settings [13, 31, 11, 23, 24, 25]. In this paper, we again perform the Strang splitting to the Vlasov equation, with which each split lower-dimensional ‘linear’ equation can be evolved by the proposed SL DG algorithm. When the relativistic Vlasov equation [18] is considered, the split 1-D equation in the form of equation (1.1) could be with general variable coefficients a(x, t). a(x, t) is reduced to a constant for the non-relativistic Vlasov equation. The SL approach has a fixed set of numerical meshes as the Eulerian approach; however in each of the time step evolution, information propagates along the characteristics in the Lagrangian fashion. Therefore, the SL method can be designed to be highly accurate as the Eulerian approach [37, 17], while being free of CFL time step restriction, owing to its Lagrangian evolution mechanism. Because of these properties, the SL approach has been very popular in the weather prediction [29] and Vlasov simulation community in the past decade [28, 2, 3, 18, 31, 4, 23]. Depending on the solution space, there have been developments of finite difference [4, 23, 24, 25], finite volume [13, 31, 11, 25] and finite element [12, 22, 30, 27, 26] SL methods in terms of algorithm design, analysis and applications. This paper focuses on the finite element SL method, which can be generally characterized by the following. 1. A finite element solution space, either continuous [12, 22, 30, 32] or discontinuous [27, 26] finite element solution space. In this paper, we consider discontinuous piecewise 2

polynomial solution space, Vhk = {v : v|Ii ∈ P k (Ii ),

Ii

is a discretized cell with length bounded by h} (1.3)

where P k (Ii ) denotes the set of polynomials of degree at most k on the cell Ii . The numerical solution at some discrete time tn is denoted by unh . 2. Exact evolution of numerical solution for a numerical time step ∆t. In this paper, the exact solution evolution operator is denoted as S∆t . 3. Projection of the exactly evolved numerical solution S∆t (unh ) back onto the finite element solution space. The most distinctive properties of finite element methods, when compared with finite volume and finite different methods [6, 28, 4, 23, 24], include their flexibility, h-p adaptivity, compactness and ability in handling complicated geometry. The use of discontinuous piecewise polynomial solution space Vhk shows great advantage over the traditional conforming finite element method in its flexibility and compactness. The proposed SL DG method is formulated based on time integrating the semi-discrete DG method proposed by Cockburn and Shu [8]. The formulation is similar in the spirit to the method proposed by Restelli et. al. in [26]. The essential difference is the exact versus numerical implementation of the solution evolution operator S∆t (unh ). The solution evolution in [26] is realized by numerical integration, which is not as accurate and introduces numerical instability when the time step is large (see Example 5.2 in Section 4). Based on SL DG, a recently developed maximum principle preserving (MPP) or positivity preserving (PP) limiter [35, 36] is applied to guarantee the MPP or PP property of the numerical solution while maintaining its high order accuracy whenever the analytical solution enjoys the MPP or PP property. The algorithm is extended to 2-D problems via Strang splitting. To guarantee the high order spatial accuracy of the scheme, we propose to work with nodes (k + 1 point values at Gaussian points per cell) rather than modes (up to k th moments per cell) in the 2-D setting to accommodate the shearing of advection velocities over a rectangular cell. We analyze the proposed SL DG method by establishing its connection with the direct and weak forms of Lagrangian Galerkin methods. In [22], direct and weak formulations of 3

the Lagrangian Galerkin method are presented for the advective and conservative form of the equations respectively. These two formulations are equivalent when the velocity field a(x, t) is divergence free, but different when ∇x · a 6= 0. We show in Section 5 that the proposed SL DG method for equation (1.1) is equivalent to the weak form of the Lagrangian Galerkin formulation using piecewise polynomial solution space Vhk . Based on such connection, we analyze the mass conservation, L2 stability and accuracy of the proposed SL DG scheme. The paper is organized as follows. Section 2 is a review of the semi-discrete DG scheme for hyperbolic equations [6]. Section 3 presents the proposed PP SL DG method for 1-D and multi-D problems. We establish the equivalence between the proposed SL DG method and the weak form of the Lagrangian Galerkin method [22] in Section 4, and analyze its stability and accuracy properties. Section 5 presents the performance of the PP SL DG scheme in a range of 1-D and 2-D test problems. Section 6 is the application to the VP system. Section 7 gives the conclusion.

2

DG scheme

In this section, we briefly review the DG methods [8, 7, 9, 6], which are a class of high order numerical methods for solving hyperbolic equations. The method approximates the weak form of the scalar hyperbolic equation (1.2) by finding uh ∈ Vhk such that for each element K, d dt

Z

Z

ˆf (uh (t, x))·nv(x)ds−

uh (t, x)v(x)dx+ K

∂K

Z f (uh (t, x))·∇x vdx = 0,

∀v ∈ Vhk (2.1)

K

int(K) ext(K) where n is the outward unit normal of the edge of K, and ˆf (uh (t, x)) = ˆf (uh , uh , n)

is defined by an approximate Riemann solver that is upwind and consistent with the flux function f (u) = a(x, t)u. Equation (2.1) is usually further discretized in time by a stable time integrator, such as the third order strong stability preserving (SSP) Runge-Kutta (RK) method [16]. The most distinctive properties of DG methods include their h-p adaptivity, compactness and ability in handling complicated geometry, compared with the finite difference WENO scheme [19]. On the other hand, the numerical time step of the DG methods is restricted by the relatively small CFL numbers according to linear stability analysis, leading ∆t to roughly |c| ∆x ≤

1 2k+1

for schemes with Vhk solution space, where |c| is the wave speed 4

associated with a given problem [10, 20]. In this paper, we also propose a scheme based on the semi-discrete DG discretization (2.1), but we propose to integrate in time in a SL fashion.

3

PP SL DG scheme

In this section, we introduce the proposed SL DG method for the 1-D advection problem (1.1) with the MPP or PP limiter to preserve the maximum principle or positivity of the numerical solution. We extend the algorithm to 2-D problems via Strang splitting.

3.1

SL DG method for 1-D advection equation

The proposed SL DG method for the equation (1.1) is designed based on equation (2.1), Z Z d uh (t)·vdx− a(x, t)uh ·vx dx+a(x, t)uh ·v|xi+ 1 −a(x, t)uh ·v|xi− 1 = 0, ∀v ∈ Vhk . (3.1) 2 2 dt Ii Ii Notice that at each interface xi+ 1 , uh should be chosen as a single valued numerical flux 2

according to upwinding, and the test function v should be chosen as v − at xi+ 1 and as v + 2

at xi− 1 in (3.1). Rather than applying the popular SSP RK methods for the equation (3.1), 2

we propose a SL integration strategy. More specifically, we integrate equation (3.1) in time from [tn , tn+1 ] and obtain, Z Ii

un+1 h

Z

unh

· vdx =

Z

tn+1

− tn

Z

· vdx + tn

Ii

Z

tn+1



Ii

a(x, t)St (unh )

. = T0 + T1 + T2 ,

a(x, t)St (unh ) · vx dxdt · v|x−

∀v ∈ Vhk

1 i+ 2

−

a(x, t)St (unh )

 · v|x+

1 i− 2

dt, (3.2)

where St (unh ) stands for the exact evolution operator starting from the initial condition unh at the time level tn . We evaluate T1 and T2 in a SL fashion. Specifically, we observe in Figure 3.1 that at a fixed spatial location at time level tn+1 , say (xi−1/2 , tn+1 ), there exists a backward characteristic line (or curve when a is a function of x and t), with the foot located on the time level tn at x?i−1/2 . Let Ωi−1/2 denote the region bounded by the three points (xi−1/2 , tn+1 ), (xi−1/2 , tn ) and (x?i−1/2 , tn ), we apply the Divergence Theorem to the integral

5

form of equation (1.1) over the region Ωi−1/2 , Z Z n n 0= (St (uh ))t + (a(x, t)St (uh ))x dxdt = Ωi−1/2

St (unh ) · nt + a(x, t)St (unh ) · nx ds. (3.3)

∂Ωi−1/2

Hence, Z

tn+1

Z

a(x, τ )Sτ (unh )(xi−1/2 , τ )dτ

tn

xi−1/2

= x?i−1/2

unh (ξ)dξ,

(3.4)

and Z

xi+ 1

2

T2 = x? 1 i+ 2

unh dx

·

) v(x− i+ 12

Z −

xi− 1 2

x? 1 i− 2

). unh dx · v(x+ i− 1 2

Similarly, T1 can be evaluated by first applying Gaussian quadrature rules, then applying the SL idea as in (3.4), Z

Z

a(x, t)St (unh )dt

T1 = Ii

= ∆xi

!

tn+1

tn

X

Z

tn

= ∆xi

ig

Z

!

tn+1

wˆig

ig

X

· vx dx

a(x, t)St (unh )dt|x=ˆxig

· vx |x=ˆxig

x ˆig

wˆig x ˆ?ig

uh (ξ, tn )dξ · vx |x=ˆxig ,

(3.5)

where xˆig and wˆig are the corresponding Gaussian quadrature nodes and weights over each of the subinterval Ii and xˆ?ig is the foot of the characteristics at tn from (ˆ xig , tn+1 ). This completes the description of the proposed SL DG scheme for the 1-D equation (1.1). The proposed SLDG scheme enjoys the advantages that (1) the scheme is of conservative form, where the flux function is evaluated exactly by equation (3.4); (2) the scheme does not have the CFL/linear stability time step restriction as the Eulerian RK DG scheme [8], leading to significant computational savings; (3) when compared with the finite difference WENO scheme, the proposed algorithm is more compact and allows for h-p adaptivity. Remark 3.1. The proposed SL DG framework is similar to the one proposed in [26]. The essential difference is that we use the Divergence Theorem to evaluate the terms T1 and T2 in equation (3.2) exactly in time via equation (3.4), while a backward characteristic tracing procedure together with Gaussian quadrature formula is used to evaluate them in [26]. The proposed SL DG scheme in this paper not only is more accurate (exact versus numerical 6

v rrr

v

v

v

v

v rrr

v

tn+1

v

v

v rrr

v

tn

       v r r r v

Ωi−1/2 v

x1/2 6 xi−5/2 xi−3/2 xi−1/2 xi+1/2 xi+3/2 x?i−1/2

xN +1/2

Figure 3.1: SL scheme approximates equation (1.1). integration) but also more robust and stable. Especially, when the time step is large (larger than the CFL restriction as in the Eulerian method), the Gaussian quadrature formula will not give an accurate approximation, as the exactly evolved solution St (unh ) is discontinuous over [tn , tn+1 ]. Moreover, numerical solution is observed to be unstable demonstrated in Figure 5.2 for Example 5.2 in Section 4.

3.2

MPP and PP limiters

In [35, 36], a MPP or PP limiter is introduced for the RK DG scheme to preserve the maximum principle or positivity, whichever is relevant for the exact solution, of the numerical solutions for hyperbolic type PDEs. The procedure of the MPP or PP limiter can be viewed as controlling the maximum, minimum, or positivity of the numerical solution (polynomials on discretized cells) by a linear rescaling around cell averages, with the assumption that the cell averages are well bounded by the maximum, minimum, or positivity enjoyed by the exact solution. In particular, for MPP, we would like to modify the numerical solution pk (x) to p˜k (x), approximating a function u(x) on a cell Ii , such that it satisfies • Accuracy: for smooth function u(x), kpk (x) − p˜k (x)k = O(∆xk+1 ). • Conservativity:

R Ii

p˜k (x)dx =

R Ii

pk (x)dx.

• MPP: p˜k (x) ∈ [m0 , M0 ], where m0 and M0 are the minimum and maximum of u(x) respectively. 7

In order to achieve the above mentioned properties, one can apply the following limiter   M0 − p¯ m0 − p¯ |, | |, 1 , (3.6) p˜k (x) = θ(pk (x) − p¯) + p¯, θ = min | 0 M − p¯ m0 − p¯ where M 0 and m0 are the maximum and the minimum of pk (x), if p¯ ∈ [m0 , M0 ]. It can be easily checked that with the application of such limiter, the conservativity and MPP properties of the numerical solution are satisfied. Furthermore, it was proved [34, 36] that such limiting process maintains the original (k + 1)th order accuracy of the approximation. This limiter can be applied to the proposed SL DG scheme. In the context of solving the Vlasov equation for example, we have m0 = 0 and M0 = 1. It is proved later in Section 4 that the cell averages of the numerical solution p¯ is the projection of the zeroth moment of the exactly evolved solution, therefore p¯ ∈ [0, 1]. If the exact solution satisfies only the positivity property, a similar PP limiter [36], 

p˜k (x) = θ(pk (x) − p¯) + p¯,

 p¯ θ = min | 0 |, 1 , m − p¯

(3.7)

can be used to enforce the positivity preserving property while maintaining the original high order accuracy. To implement the MPP (3.6) or PP (3.7) limiter, the maximum and/or minimum of the numerical solution pk (x) (m0 and/or M 0 ) are needed. In our numerical tests in Section 5 and 6, we use up to P 3 polynomials, the maximum and the minimum of which are easily found by locating the zeros of their derivatives. We can also modify the MPP or PP limiter, as was done for the RK DG method in [35, 36], so that only the maximum and minimum of the solution at certain Gaussian points are used. We should however notice that the integrand unh (ξ) in (3.4) over the integral region (x?i−1/2 , xi−1/2 ) may be a piecewise polynomial rather than a single polynomial, hence the choice of these Gaussian points should be over each sub-interval within (x?i−1/2 , xi−1/2 ) on which unh (ξ) is a polynomial, so that the quadrature computes the integral in (3.4) exactly. This strategy is useful for very high order polynomials where the exact maximum and minimum of pk (x) are no longer easy to find.

3.3

2-D implementation

To extend to 2-D problems, we adapt the idea of Strang splitting which is O(∆t2 ) accurate. Because of the splitting, our 2-D numerical mesh is the tensor product of two 1-D meshes, 8

u

u

u

u

u

u

u

u

u

u

u

u

u

u

u

u

-

-

Figure 3.2: k = 3. therefore a rectangular mesh. It is important to take the shearing of advection velocities over a cell into account, in order to design a scheme that has high order spatial accuracy. Because of this, we consider the solution space as (k + 1)2 point values at Gaussian nodes per cell. Take the rotational problem ut − yux + xuy = 0 for example, the 2D algorithm based on the 1D SL DG algorithm with solution space Vhk via splitting is the following. 1. In each of the rectangular cell, locate (k + 1) Gaussian quadrature nodes in both x and y directions as (xig , yjg ), ig, jg = 1, · · · k + 1. For example, see Figure 3.2 for the case of k = 3. 2. Evolve ut − yux + xuy = 0 by Strang splitting. For the split equation, e.g. ut − yux = 0, evolve the 1-D problems at different yjg ’s with different advection velocities, see Figure 3.2. For each yjg , the (k + 1) point values are mapped to a P k polynomial per cell (the unique interpolation polynomial of degree up to k); then the 1-D problem is evolved by the proposed SL DG scheme; finally the evolved P k polynomial is mapped back to the (k + 1) point values at Gaussian nodes to update the solution.

4

Theoretical properties

In this section, we first review the direct and weak forms of the Lagrangian Galerkin method [22]. Then we establish the connection between the SL DG and the weak form of the Lagrangian Galerkin method, assuming the integrals in equation (3.2) are evaluated exactly. Based on such connection, we analyze the L1 conservation, L2 stability and error estimate for the proposed SL DG scheme. We remark that although the SL DG scheme is based on the 9

semi-discrete form of the original DG scheme, for which the stability and accuracy property has been well established [6], the corresponding results for the SL DG scheme can not be concluded directly. The essential difference, for example in the stability analysis for the semi-discrete version of the DG scheme, is that the test function v in equation (3.1) is taken as the numerical solution uh (x, t) which is time dependent. However, for the SL DG scheme, the test function v is time independent, see equation (3.2). Toward the end of this section, we drop the assumption of exact integration, but consider using the Gaussian quadrature in evaluating T1 in (3.2). We show that the proposed SL DG scheme is reduced to the exact shifting of moments when CF L = 1, assuming a uniform mesh and constant advection speed. This property can be used to show that, when CF L is large, only local information around the domain of dependence at tn is actually involved in updating the solution, despite the appearance of ‘global information involvement’ in the scheme formulation (3.2).

4.1

Direct and weak forms of the Lagrangian Galerkin method

Consider the scalar linear advection equation in advective and conservative form respectively, ut + a · ∇x u = 0,

(4.1)

ut + ∇x · (au) = 0.

(4.2)

When the velocity field a is divergence free, i.e. ∇ · a = 0, the above two equations are equivalent. We define the characteristic paths/trajectories X(x, s; t) by Z t a(X(x, s; τ ), τ )dτ. X(x, s; t) = x +

(4.3)

s

Especially, X(y, tn+1 ; tn ) and X(x, tn ; tn+1 ) will be denoted by x and y respectively. In the following, (·, ·) denotes the inner product of two functions. The direct Lagrangian Galerkin in terms of formulation in evolving the advective form of the equation (4.1) from unh to un+1 h finite element basis function Φj is (un+1 h , Φj )

Z =

unh (x)Φj (y)dy,

∀j = 1, · · · , Nb ,

(4.4)

where Nb is the number of basis functions. The weak Lagrangian Galerkin formulation for the conservative form of the equation (4.2), on the other hand, is Z n+1 (uh , Φj ) = unh (x)Φj (y)dx, ∀j = 1, · · · , Nb . 10

(4.5)

Remark 4.1. The essential difference between the two formulations above is that in the direct formulation (4.4), the information of the numerical solutions unh is being tracked forward along characteristics; whereas in the weak formulation, the information for the test function Φj is being tracked backward along characteristics. It was commented in [22] that these two formulations are equivalent when ∇ · a = 0, but are different otherwise and are designed for the advective and conservative forms of the equation (4.1) and (4.2) respectively. Remark 4.2. Originally developed for the continuous finite element solution space [22], the above direct and weak formulations can be generalized to DG solution space Vhk , where Φj is a polynomial basis function on a discretized cell.

4.2

Connection between the SL DG method and the weak form of the Lagrangian Galerkin method

In this subsection, we establish the connection between the proposed SL DG method and the Lagrangian Galerkin formulations reviewed above in the 1-D case. Especially, assuming that the integrals in equation (3.2) are evaluated exactly, we show that the SL DG method (3.2) is equivalent to both the direct and the weak formulations (4.4) and (4.5), when a is a constant, with the DG solution space Vhk . We then show that for general variable coefficient a, the proposed SL DG method is only equivalent to the weak formulation (4.5). Proposition 4.3. The SL DG algorithm (3.2) is equivalent to the direct formulation (4.4) when a is constant. Proof. It is sufficient to show the equivalency for v = Φi,l , where Φi,l ∈ P k (Ii ), ∀i, l. Performing integration by parts on equation (3.2) gives, (un+1 h , Φi,l )

=

(unh , Φi,l )

Z

tn+1

Z

− tn

= (unh , Φi,l ) + = (unh , Φi,l ) +

Z

tn+1

Ii

Z

tn

Z Ii

11

(Sτ (unh ))τ Φi,l dxdτ

Ii

(S∆t (unh ) − unh )Φi,l dx

= (S∆t (unh ), Φi,l ) Z = unh (x)Φi,l (y)dy. Ii

(aSτ (unh ))x Φi,l dxdτ

Proposition 4.4. The SL DG algorithm (3.2) is equivalent to the weak formulation (4.5) when a is constant. Proof. Since the direct and weak formulations (4.4) and (4.5) are equivalent when a is constant, it is straightforward to make the claim. However, in the following, we will prove it in a way that can be generalized to a variable coefficient a. We consider equation (3.2) with v ∈ P k (Ii ). Without loss of generality, we assume CF L < 1 and a > 0. When CF L ≥ 1, the scheme can be decoupled to the situations with CF L < 1. We introduce the following notations, . n n Xi± 1 (t) = X((xi± 1 , t ); t),

. n+1 n+1 Xi± ); t), 1 (t) = X((xi± 1 , t

2

2

t ∈ [tn , tn+1 ]

2

2

with the notation for characteristics (see equation (4.3)) and let . X?,i± 1 = X((xi± 1 , tn+1 ); tn ).

. ? n n+1 Xi± ), 1 = X((xi± 1 , t ); t 2

2

2

2

Consider the term T1 in equation (3.2), Z Z Z  Z . n T1 = aSτ (uh )·vx dxdτ = + + aSτ (unh )·vx dxdτ = T11 +T12 +T13 (4.6) I1

Ii ×[tn ,tn+1 ]

where Z

Z

tn+1

Z

Xn

= tn

I1

i− 1 2

(τ )

Z ,

I2

Z

tn+1

I3

Z

X n+1 1 (τ ) i+ 2

= I2

xi− 1

Z ,

2

tn+1

=

X n 1 (τ ) i− 2

tn

Z

I3

tn

Z

xi+ 1 2

X n+1 (τ ) i+ 1 2

are the three regions in the rectangle Ii × [tn , tn+1 ] divided by the two characteristic lines n+1 n Xi− (t), ∀t ∈ [tn , tn+1 ]. First we claim that 1 (t) and X i+ 1 2

2

Z

X?,i+ 1 2

T12 =

unh (x)v(y)dx

Z −

xi− 1

X?,i+ 1

2

unh (x)v(x)dx.

(4.7)

xi− 1

2

2

To see this, we perform a change of variables from (x, t) in the Cartesian coordinates to (ξ, η) in the characteristic coordinates, where x = ξ + a(η − tn ),

12

t = η.

(4.8)

Then Z

tn+1

Z

X?,i+ 1 2

T12 = tn

aunh (ξ)vξ dξdη

(4.9)

unh (ξ)vη dξdη

(4.10)

xi− 1

2

Z

tn+1

Z

X?,i+ 1 2

= tn

xi− 1

2

Z

X?,i+ 1

2

=

unh (ξ)

!

tn+1

Z

vη dη dξ tn

xi− 1 2

Z

X?,i+ 1 2

=

unh (ξ)(v(ξ + a∆t) − v(ξ))dξ

xi− 1 2

Z

X?,i+ 1

2

=

unh (x)v(y)dx

2

∂v ∂x ∂x ∂η

where the second equality above is due to vη = claim that 2

T11 =

unh (x)v(x)dx,

xi− 1

2

xi− 1

2

−

xi− 1

Z

X?,i+ 1

Z

unh (x)

∂v ∂v ∂x = a ∂x = a ∂x = avξ . Secondly, we ∂ξ

tn+1

Z

aSτ (unh )dτ · v|x+ 1 .

· v(y)dx −

(4.11)

i− 2

tn

X?,i− 1 2

To see this, we start with Z

tn+1

T11 =

 Z 

tn

Xn

i− 1 2

(τ )

Z −



xi− 1 2

X n+1 (τ ) i− 1 2

X n+1 1 (τ ) i− 2

 a(Sτ (unh ))(x)vx dxdτ.

Similar to (4.7), the first integral above gives Z

Xi− 1

2

unh (x)(v(y)

Xi− 1

Z

2

− v(x))dx =

X?,i− 1

unh (x)v(y)dx,

(4.12)

X?,i− 1

2

2

and the second integral gives Z

xi− 1

2

X?,i− 1

Z

T ((ξ,tn );xi− 1 ) 2

tn

unh (ξ)avξ dξdη =

Z

xi− 1 2

unh (ξ)

X?,i− 1

2

Z

T ((ξ,tn );xi− 1 ) 2

! vη dη dξ

tn

2

Z

xi− 1

2

=

unh (ξ)dξ · v|x+

X?,i− 1

1 i− 2

2

Z

tn+1

= tn

aSτ (unh )dτ · v|x+ 1 ,

(4.13)

i− 2

where T ((ξ, tn ); xi− 1 ) is the time taken for the characteristics to propagate from (ξ, tn ) to 2

xi− 1 , with the same change of variables as defined in (4.8). Combining equations (4.12) and 2

13

(4.13) gives the claim (4.11). With a similar analysis, we have Z x 1 Z tn+1 i+ 2 n T13 = − uh (x) · v(x)dx + aSτ (unh )dτ · v − |xi+ 1 .

(4.14)

2

tn

X?,i+ 1

2

Plugging equations (4.11), (4.7) and (4.14) into (4.6), and into the SL DG scheme (3.2) gives Z n+1 (4.15) (uh , v) = unh (x)v(y)dx, which is the weak formulation (4.5). Next, we generalize the equivalence established in Prop. 4.4 to equation (1.1) with a general variable coefficient a(x, t). Proposition 4.5. The SL DG algorithm (3.2) is equivalent to the weak formulation (4.5) for general variable coefficient a(x, t). Proof. The proof follows similarly as that for Prop. 4.4. We divide the rectangle Ii ×[tn , tn+1 ] by the two characteristics curves, analyze the T1 term as three separate terms T11 , T12 and T13 , and claim equations (4.11), (4.7) and (4.14). The essential difference is that the characteristics curves, therefore the coordinates transform between (x, t) and (ξ, η) are now defined by, dx(ξ, η) = a(x(ξ, η), η), dη

x(tn ) = ξ;

t = η.

(4.16)

Take the analysis for T12 for example. In order to prove equation (4.7), we need to generalize equation (4.9) and (4.10) for the case of variable coefficient as following,   Z η Z tn+1 Z X 1 ?,i+ 2 T12 = a(x(ξ, η), η) exp − ax (x(ξ, τ ), τ )dτ unh (ξ)vξ dξdη tn

tn

Xi− 1 2

Z

tn+1

Z

X?,i+ 1 2

= tn

unh (ξ)vη dξdη,

Xi− 1 2

where the first equality is due to the fact that u(x, t) exp

R

t tn

ax (x, τ )dτ



is constant along

the characteristic curves. To show the second equality, it is sufficient to show  Z η  ax (x, τ )dτ xξ = xη . a(x, η) exp −

(4.17)

tn

To see this, we know xη = a(x, η) and xξ = 1 + dxξ = ax xξ dη

⇒

Rη

ax xξ dη from equation (4.16). Therefore, Z η  xξ = exp ax dη , tn

tn

14

hence equation (4.17). Similar modification can be applied to prove equations (4.11) and (4.14), therefore the weak formulation (4.5). Remark 4.6. The equivalence of the SL DG scheme (3.2) and weak formulation (4.5) is established based on the assumption that the integrals in both formulations are evaluated exactly. When the integrals are approximated by numerical quadratures, the two formulations are different. Remark 4.7. When ∇x · a 6= 0, the proposed SL DG is equivalent to the weak formulation (4.5), but not the direct formulation (4.4).

4.3

Stability and accuracy analysis

The stability and accuracy analysis for linear advection equation with constant coefficient a are based on the weak formulation (4.5), which were presented in [22] for conforming finite elements. Proposition 4.8. (L1 (mass) conservation) The PP SL DG scheme is mass conservative, if periodic boundary condition is imposed Z Z n+1 uh dx = unh dx. I

(4.18)

I

Assuming that u(x, t = 0) ≥ 0, when the PP limiter is applied, the scheme is conservative in the L1 norm, n kun+1 h k1 = kuh k1 .

(4.19)

Proof. Let the test function v = 1 in equation (3.2) and sum the equation up over i gives the mass conservation equation (4.18), due to the cancellation of flux terms between neighboring cells. In fact, equation (3.2) with v = 1 gives local conservation (changes of the mass in a local cell come from the flux difference at the cell boundaries with locally defined fluxes), which is a stronger property than the global mass conservation (4.18). The PP limiter preserves such conservation [36]. Moreover, the PP limiter preserves the positivity of the solution, i.e. unh ≥ 0 when u(x, t = 0) ≥ 0, hence equation (4.19). Notice that this proposition holds for equation (1.1) with variable coefficient a(x, t) as well.

15

Proposition 4.9. (L2 stability) Consider the SL DG scheme (3.2), or its equivalent weak form (4.5), with the solution space Vhk for the linear advection equation (1.1) with constant coefficient a and periodic boundary condition. Assume that the integrals in equation (3.2) are being exactly evaluated, then n kun+1 h k2 ≤ kuh k2 .

(4.20)

Proof. We prove the L2 stability via the weak formulation (4.5). Especially, let the test function Φj = un+1 h , n+1 n+1 n (un+1 h , uh ) = (uh , S−∆t (uh ))

≤ kunh k2 kS−∆t un+1 h k2 ≤ kunh k2 kun+1 h k2 therefore L2 stability (4.20). Proposition 4.10. (error estimate) The error for the proposed SL DG scheme (3.2) or its equivalent form (4.5) with solution space Vhk for the linear advection equation (1.1) with constant coefficient a and periodic boundary condition satisfies kunh − u(x, tn )k2 ≤ CN ∆xk+1 , where N is the number of time steps taken. Especially when ∆t = O(∆x) and T = N ∆t is a constant, the error analysis indicates k th order convergence, which is suboptimal. Proof. Let en = u(x, tn ) − unh ,

η n = u(x, tn ) − Πu(x, tn ),

ξ n = en − η n ,

where Πu is defined to be the L2 projection of u into the finite element space. Since (u(x, tn+1 ), Φj ) = (u(x, tn ), S−∆t Φj ),

n (un+1 h , Φj ) = (uh , S−∆t Φj ),

subtract one equation from another gives, (en+1 , Φj ) = (en , S−∆t Φj ). 16

Let Φj = ξ n+1 , then (ξ n+1 , ξ n+1 ) = (ξ n + η n , S−∆t ξ n+1 ) ≤ (kξ n k2 + kη n k2 )kξ n+1 k2 . Therefore kξ n+1 k2 ≤ (kξ n k2 + kη n k2 ) ≤ · · · ≤

X

kη n k2 = CN ∆xk+1

n

where C depends only on the exact solution u and its derivatives.

4.4

Local/global involvement of information when CF L > 1

In the previous subsection, we assume that the integrals in the SL DG scheme (3.2) are evaluated exactly, based on which we establish the equivalence of the proposed SL DG scheme with the weak formulation (4.5). Although information appears to be involved globally in the SL DG formulation (3.2) for large time step evolution, only the information around the foot of characteristics is really involved, see equation (4.5). However, exact evaluation of integrals (T1 in (3.2)) is in general very difficult for variable coefficient a(x, t). In this subsection, we show the local involvement of information when T1 is evaluated by Gaussian quadrature (see (4.6)) for large time step evolution, assuming the numerical mesh is uniform. Proposition 4.11. When CF L = 1, the proposed SL DG scheme for linear advection equation ut + ux = 0 is reduced to exact shifting of moments, when the numerical mesh is uniform for the solution space Vhk , ∀k. Proof. We first introduce the following notations. A discretized cell is Ii = [xi− 1 , xi+ 1 ] with 2

2

midpoint xi and interval length ∆xi . The orthogonal basis functions on Ii are {Φl (ξ)}kl=0 R1 . i with −2 1 Φl (ξ)Φm (ξ)dξ = δlm where ξ(x) = x−x ∈ [− 21 , 21 ] for the Vhk space. The numerical ∆xi 2

solution on Ii written as ui (x, t) =

k X

uˆi,l (t)Φl (ξ(x)),

x ∈ Ii

l=0

satisfies equation (3.2) for any v = Φm , m = 0, · · · k, i.e. Z Z x n+1 n uˆi,m = uˆi,m + u(ξ, tn )dξ · Φ0m (ξ(x))dx Ii

Z

xi+1/2

x?

n

Z

xi−1/2

u(ξ, t )dξ · Φm (1/2) −

− x?i+1/2

! n

u(ξ, t )dξ · Φm (−1/2) ,

m = 0, · · · k.

x?i−1/2

(4.21) 17

Especially, when CF L = 1 (∆t = ∆x), equation (4.21) becomes, Z Z x n+1 n uˆi,m = uˆi,m + un (η)dη · Φ0m (ξ(x))dx − ∆x · u¯ni−1 · Φm | 1 + ∆x · u¯ni · Φm |− 1 Ii (4.6)

uˆni,m + ∆x

=

2

x−∆x

X

Z

ξ −1 (ˆ x

ig )

un (η)dη · Φ0m (ˆ xig ) − ∆x · u¯ni−1 Φm | 1 + ∆x · u¯ni Φm |− 1

wˆig ξ −1 (ˆ x

ig

2

2

ig )−∆x

∀m = 0, · · · k,

2

(4.22)

where ξ −1 (ˆ x) = xi + xˆ∆xi . Consider the second term in the equation above, ∆x

X

Z wˆig

ξ −1 (ˆ xig )−∆x

ig

= ∆x

X

ξ −1 (ˆ xig )

wˆig

= ∆x

X

xi+ 1

Z Z ( − Z

wˆig

(ˆ u0,i −

Z

2

ξ −1 (ˆ x

Ii

ig

u(η, tn )dη · Φ0m (ˆ xig ) +

2

ξ −1 (ˆ x

ig )

!

xi− 1

)u(η, tn )dη

· Φ0m (ˆ xig )

ig )−∆x

!

xi+ 1

2

(u(η, tn ) − u(η − ∆x, tn ))dη)

ξ −1 (ˆ xig )

ig

Z Z

1 2

= ∆x · uˆ0,i Φm (x)|− 1 − 2

1

= ∆x · uˆ0,i Φm (x)|−2 1 2

xi+ 1 2

n

n

· Φ0m (ˆ xig ) 

(u(η, t ) − u(η − ∆x, t ))dη dΦm (x) Z + ∆x · (ˆ u0,i − uˆ0,i−1 )Φm (x)|− 1 − (u(x, tn ) − u(x − ∆x, tn ))Φm (x)dx Ii

x

2

= ∆x · uˆ0,i Φm (x)| 1 − ∆x · uˆ0,i−1 Φm (x)|− 1 − 2

2

(ˆ uni,m

Ii

− uˆni−1,m ),

where the third equation from the end is due to the exactness of Gaussian integration for polynomials and the second equation from the end is due to integration by parts. Plugging this into equation (4.22) gives uˆn+1 ˆni−1,m . i,m = u Corollary 4.12. Assume the numerical mesh is uniform, when CF L > 1, only local information around the domain of dependence is involved in updating the solution for ut +ux = 0, despite the appearance of ‘global involvement of the information’ in the proposed SL DG scheme for any Vhk . Proof. The proof is based on Prop. 4.11 and the application of the fact Z tn+1 Z tn+1 ! J Z tn +j∆x X g(x)dx = + g(x)dx, ∀g(x), where J = bCF Lc tn

j=1

tn +(j−1)∆x

tn +J∆x

to equation (3.2). We omit the details.

18

Remark 4.13. Unfortunately, the property of the ‘local’ involvement of information around the domain of dependence as described in Cor. 4.12 does not apply to the case of non-uniform meshes.

5 5.1

Numerical tests Basic 1-D and 2-D tests

Example 5.1. (one dimensional linear advection) ut + ux = 0,

x ∈ [0, 2π].

(5.1)

For this example the exact solution satisfies the maximum principle and hence the MPP limiter is relevant. The proposed MPP SL DG methods are used to solve equation (5.1) with smooth initial data u(x, 0) = sin(x) and exact solution u(x, t) = sin(x − t). The L2 error and the corresponding order of convergence for schemes with uniform and nonuniform meshes and different solution spaces Vhk are summarized in Tables 5.1 and 5.2. We also apply the scheme to advect a rectangular wave, see Figure 5.1. As the MPP limiter is applied here to enforce a maximum principle, the numerical solutions are observed to be bounded by the lower and upper bounds of the initial condition. However, we observe that oscillations are not completely removed near the discontinuities. This is consistent with the observation in [36], that the MPP limiter should be used to augment rather than to replace high order preserving non-oscillatory limiters such as the total variation bounded (TVB), essentially non-oscillatory (ENO) or weighted ENO (WENO) type limiters, when the complete removal of numerical oscillations is desired. Example 5.2. (one dimensional equation with variable coefficient) ut + (sin(x)u)x = 0 x ∈ [0, 2π].

(5.2)

The initial condition is u(x, 0) = 1 and the boundary condition is periodic. The exact solution is
sin 2 tan−1 (e−T tan( x2 )) u(x, t) = . sin(x) For this example the exact solution satisfies the positivity preserving property but not the maximum principle, hence the PP limiter is relevant. 19

Table 5.1: MPP SL DG methods with different solution spaces Vhk for solving (5.1) with u(x, t = 0) = sin(x) at T = 20. CF L = 2.2. The numerical meshes are uniform. L2 error mesh 20 40 80 160 320

k=0 error order 3.96E-1 – 2.08E-1 0.92 1.06E-1 0.96 5.39E-2 0.98 2.71E-3 0.99

k=1 error order 1.61E-2 – 2.70E-3 2.11 6.48E-4 2.06 1.54E-4 2.07 3.73E-5 2.05

k=2 k=3 error order error 2.10E-4 – 4.28E-6 2.60E-5 3.02 2.50E-7 3.23E-6 3.00 1.49E-8 4.04E-7 3.00 9.09E-10 5.05E-8 3.00 5.59E-11

order – 4.09 4.07 4.04 4.02

Table 5.2: MPP SL DG methods with different Vhk for solving (5.1) with u(x, t = 0) = sin(x) at T = 20. CF L = 2.2. The numerical meshes are non-uniform, with a 10% perturbation from the uniform meshes. L2 error mesh 20 40 80 160 320

k=0 error order 1.14E-0 – 8.14E-1 0.49 5.00E-1 0.70 2.79E-1 0.84 1.47E-2 0.92

k=1 error order 1.81E-2 – 4.17E-3 2.11 9.86E-4 2.08 2.30E-4 2.09 5.54E-5 2.06

k=2 k=3 error order error 5.53E-4 – 2.91E-5 6.77E-5 3.03 1.79E-6 8.37E-6 3.01 1.06E-7 1.04E-6 3.00 6.54E-9 1.30E-7 3.00 3.93E-10

order – 4.02 4.07 4.02 4.05

The proposed PP SL DG methods are applied to solve equation (5.2). Tables 5.3 gives the L2 error and the corresponding order of convergence for schemes with non-uniform meshes and different solution spaces Vhk . Expected order of accuracy is observed. We have also implemented the SL DG scheme with P 0 solution space proposed in [26] to this example. When CF L ≤ 1, the scheme works well, as it is reduced to the classical first order monotone finite volume scheme. However, when CF L > 1, the numerical solution is observed to be unstable, even when the exact solution is still smooth, see Figure 5.2. This can be explained by the fact that the solution at some cell boundaries over a time step [tn , tn+1 ] has discontinuities when CF L > 1. Numerical Gaussian quadrature rule does not provide an R tn+1 accurate approximation to tn a(x)u(xi−1/2 , τ )dτ in equation (3.4). Example 5.3. (two dimensional linear transport) ut + ux + uy = 0,

x ∈ [0, 2π], 20

y ∈ [0, 2π].

(5.3)

Figure 5.1: The numerical solution of MPP SL DG method for linear advection of rectangular wave with CF L = 2.2 and final integration time 20. The numerical mesh is uniform (left) and non-uniform with a 10% perturbation from the uniform mesh (right) with 90 cells. Table 5.3: PP SL DG scheme with different Vhk for (5.2) with u(x, t = 0) = 1 at T = 1 with CF L = 2.2. The numerical meshes are non-uniform based on a 10% perturbation of uniform meshes. L2 error mesh 20 40 80 160 320

k=0 error order 3.65E-1 – 1.77E-1 1.04 9.50E-2 0.90 4.77E-2 0.99 2.44E-2 0.97

k=1 error order 5.25E-2 – 1.01E-2 2.38 2.76E-3 1.87 6.98E-4 1.98 1.84E-4 1.92

k=2 error order 4.46E-3 – 7.26E-4 2.62 9.71E-5 2.90 1.18E-5 3.04 1.52E-6 2.95

k=3 error order 1.31E-3 – 5.17E-5 4.66 3.32E-6 3.96 2.07E-7 4.00 1.42E-8 3.87

For this example the exact solution satisfies the maximum principle and hence the MPP limiter is relevant. The equation is being split into two one-dimensional equations, each of which is evolved by the proposed MPP SL DG method. For 2-D linear transport equations, the SL method is essentially a shifting procedure. Since the x-shifting and y-shifting operators commute, there is no dimensional splitting error in time and the spatial error is the dominant error. Table 5.4 gives the L2 error and the corresponding order of convergence when the MPP SL DG schemes are applied to equation (5.3) with smooth solution u(x, y, t) = sin(x + y − 2t). In all of our 2-D simulations including VP simulations in the

21

Figure 5.2: Plots of the numerical solutions of the P 0 SL DG method proposed in [26] for equation (5.2) with CF L = 2.2 and final integration time 1. The numerical mesh is 80 (left) and 200 (right). next section,  CF L = ∆t

cx cy + min ∆x min ∆y

 ,

where cx and cy are the maximum wave propagation speed in the x and y directions respectively. Table 5.4: Order of accuracy of L2 error for solving (5.3) with u(x, y, t = 0) = sin(x + y) at T = 10. CF L = 2.2. The numerical meshes are non-uniform based on a 10% perturbation of uniform meshes. L2 error mesh 20×20 40×40 60×60 80×80 100×100 120×120 140×140 160×160

k=0 error order 5.47E-1 – 3.52E-1 0.63 2.53E-1 0.81 1.97E-1 0.87 1.62E-1 0.86 1.37E-1 0.92 1.19E-1 0.91 1.05E-1 0.95

k=1 error order 3.02E-2 – 8.02E-3 1.91 4.04E-3 1.69 2.31E-3 1.95 1.61E-3 1.62 1.11E-3 2.03 8.54E-4 1.71 6.67E-4 1.85

k=2 error order 1.18E-3 – 2.35E-4 2.33 8.99E-5 2.37 4.44E-5 2.46 2.52E-5 2.53 1.57E-5 2.61 1.04E-5 2.67 7.23E-6 2.71

k=3 error order 3.72E-5 – 3.06E-6 3.60 7.02E-7 3.63 2.08E-7 4.23 1.02E-7 3.19 5.48E-8 3.40 3.06E-8 3.76 1.95E-8 3.39

Example 5.4. (rigid body rotation) ut − (yu)x + (xu)y = 0,

x ∈ [−π, π], 22

y ∈ [−π, π].

(5.4)

Figure 5.3: Plots of the initial profile. For this example the exact solution satisfies the maximum principle and hence the MPP limiter is relevant. The equation is being Strang split into two one-dimensional equations, each of which is evolved by the proposed MPP SL DG method. The initial condition we used is plotted in Figure 5.3 in mesh and in contour. It includes a slotted disk, a cone as well as a smooth hump, similar to the one in [21] for comparison purpose. Numerical solutions after a full revolution from schemes with different solution space Vhk are plotted in Figure 5.4 by mesh and contour. Different numerical meshes are chosen for different Vhk , so that the overall degrees of freedom are the same. To look into the details of different numerical solutions and compare their performances, we plot slides of numerical solutions at different locations in Figure 5.5. The numerical solutions are observed to be bounded by the lower and upper bounds of the initial solution, due to the MPP limiter which enforces the maximum principle in this case. Sharper resolutions are observed for schemes with higher degree of polynomial in the solution space. Again, the MPP limiter does not completely remove the oscillations generated by the DG scheme. Example 5.5. (Swirling deformation flow) We consider solving x y ut − (cos2 ( ) sin(y)g(t)u)x + (sin(x) cos2 ( )g(t)u)y = 0, 2 2

x ∈ [−π, π],

y ∈ [−π, π], (5.5)

where g(t) = cos(πt/T )π and T = 1.5, the same example as in [21]. The initial condition is set the same as in Example 5.4, see Figure 5.3. The equation is being Strang split into two 23

one-dimensional equations, each of which is evolved by the proposed SL DG method. We apply the PP limiter here because the split 1-D equations are variable coefficient equations which do not have the MPP property. We set CF L = 0.2, so that the splitting error in time will not dominate. We numerically integrate the solution up to time 0.75, when the initial profile is greatly deformed, and to time 1.5, when the initial profile should be recovered. Figure 5.6 shows numerical results (in both mesh and contour plots) from the scheme at time 1.5 when the initial profile should be recovered. Figure 5.7 demonstrates the numerical results by slides benchmarked with the exact solution. The numerical solutions with high degree polynomial in the solution space are observed to better resolve the solution. Figure 5.8 shows the mesh and contour plots of the numerical solution of the scheme with P 1 solution space at integration time 0.75, when the solution is quite deformed. Oscillations are observed around the discontinuities. However, the numerical solutions are observed and verified to be non-negative with the PP limiter. We notice that the exact solution of the PDE (5.5) does satisfy the MPP property because the variable coefficients are divergence-free. If we use the MPP limiter in the SL DG method, we do get better solution as shown in Figure 5.9, 5.10, 5.11, in which the oscillations are much less than those in Figure 5.6, 5.7, 5.8. However, since in each split step the one dimensional variable coefficient PDE does not enjoy the MPP property, enforcing the MPP limiter in each split step might affect the high order accuracy of the SL DG method.

6

Vlasov Simulations

In this section, we apply the proposed MPP SL DG method to the VP system ∂f + v · ∇x f + E(t, x) · ∇v f = 0, ∂t E(t, x) = −∇x φ(t, x),

−4x φ(t, x) = ρ(t, x),

(6.1) (6.2)

which describes the collisionless plasma. In equations (6.1) - (6.2), x and v are the coordinates in the phase space (x, v) ∈ R3 × R3 , E is the electric field, φ is the self-consistent electrostatic potential and f (t, x, v) is the probability distribution function which describes the probability of finding a particle with velocity v at position x and at time t. The probability distribution function couples to the long range fields via the charge density, 24

ρ(t, x) =

R R3

f (t, x, v)dv − 1, where we take the limit of uniformly distributed infinitely

massive ions in the background. Equations (6.1) and (6.2) have been nondimensionalized so that all physical constants are one. The VP system describes the following physical process: electrons are moving with velocities v in a constant ion background; at the same time, their velocities v are being accelerated or decelerated by self-induced electric field E determined by the Poisson equation. The Strang splitting SL method for the VP system was originally proposed in [5], and soon gained wide popularity [13, 31, 11]. The Strang splitting reduces the high dimensional nonlinear Vlasov equation into lower-dimensional linear advection equations, allowing the direct application of the MPP SL DG method. The splitting decouples the simultaneous ‘particle moving’ and ‘velocity accelerating/decelerating’ processes as two separate processes: first particles are moving with constant speed v, then the particles do not move while their velocities are accelerated/decelerated by the self-induced electric field determined by Poisson’s equation. The time splitting form of equation (6.1) is, ∂f + v · ∇x f = 0 , ∂t

(6.3)

∂f + E(t, x) · ∇v f = 0 . ∂t

(6.4)

The split form of equation (6.1) can be made second order accurate in time by solving equation (6.3) for a half time step, then solving equation (6.4) for a full time step, followed by solving equation (6.3) for a second half time step. The observation that both equation (6.3) and equation (6.4) are linear hyperbolic equations makes the Strang splitting SL approach very popular in the plasma simulation community [28, 15, 33, 4, 11], as the SL scheme for the split linear equation is usually simple, effective and free of CFL restriction. The MPP SL DG method offers an alternative approach in the family of Strang split SL approaches for the VP system, since the exact solution satisfies the MPP property. For the relativistic Vlasov equation, for which only the PP property is enjoyed by the exact solution, the SL DG method with PP limiter will be used instead. Below, we briefly recall some classical preservation results in the VP system. We hope that our numerical solutions can preserve these analytically conserved quantities as much as possible. 25

1. Preservation of the Lp norm, for 1 ≤ p < ∞. Z Z d f (x, v, t)p dxdv = 0 dt v x 2. Preservation of the entropy Z Z d f (x, v, t) log(f (x, v, t))dxdv = 0 dt v x 3. Preservation of the energy Z Z  Z d 2 2 f (x, v, t)v dxdv + E (x, t)dx = 0 dt v x x

(6.5)

(6.6)

(6.7)

In this section, the proposed scheme is tested in the classical problems in plasma physics, such as Landau damping and two stream instabilities. The performance of the schemes will be demonstrated by the solution profiles, as well as by tracking the time evolution of theoretically preserved quantities (equation (6.5) - (6.7)) in the discrete sense. Periodic boundary conditions are imposed in the x-direction and zero boundary conditions are imposed in the v-direction for all of our test problems. A local DG method [1] is applied to solve the Poisson’s equation. In this paper, we only consider the 1-D Vlasov equation (6.1) with (x, v) ∈ R × R. The extension to higher dimensions in x and v can be performed by dimensional splitting without splitting error, see Example 5.3. Unless otherwise specified, our numerical simulation parameters are Nx = 64, Nv = 128 and CF L = 3, where Nx is the number of grid points along the x axis, Nv is the number of grid points along the v axis for different solution spaces Vhk . Higher degree polynomial spaces have more degrees of freedom per cell, therefore better resolution of the solution is expected. Example 6.1. (Weak Landau damping) Consider the weak Landau damping for the VP system, with the initial condition,  2 1 v f (x, v, t = 0) = √ (1 + α cos(kx)) exp − , 2 2π

(6.8)

with α = 0.01, k = 0.5 and vmax = 5, where vmax is the maximum velocity on the phase space mesh. The system is evolved by Strang split MPP SL DG method. The numerical solutions are verified to be strictly within [0, 1], because of the MPP limiter that we have 26

applied which enforces the MPP property in this case. The time evolution of the L2 norm of the electric field is plotted in the first row of Figure 6.1. The correct damping of the electric field is observed in the plots, benchmarked with the theoretical value γ = 0.1533 [14] (the solid line in the same plots) up to some time. The MPP SL DG scheme with higher degree of polynomial space performs better in matching with the theoretical damping rate in the long run, which is due to the better resolution of the MPP SL DG scheme with the same mesh. The time evolution of the L1 , L2 solution norms, energy, entropy in the discrete sense are demonstrated in the mid and bottom plots in Figure 6.1. The better performance of the scheme with higher degree polynomial space is observed. Despite this, we remark that in the weak Landau damping case, the relevant physical norms are preserved pretty well for all of the schemes we tested (see the magnitude variance in all of the y-axis) in Figure 6.1. Example 6.2. (Strong Landau damping) The next example we consider is the case of strong Landau damping. We simulate the VP system with the initial condition in equation (6.8) with α = 0.5, k = 0.5 and vmax = 5. In the first row of Figure 6.2, the time evolution of the L2 norm of the electric field is plotted for different Vhk , which is observed to be quite consistent with each other. No significant difference among the three curves is observed. The discrete L1 norm, L2 norm, kinetic energy and entropy for the MPP SL DG schemes with different Vhk spaces are also plotted in Figure 6.2. It is observed that schemes with higher degree of polynomial space Vhk in general perform better in preserving those physical norms, because of the better resolution they have per cell. Since the numerical solutions we obtain are consistent with Figure 4.9 in [23], we skip the detailed demonstration of them to save space. Example 6.3. (Two stream instability [14]) Consider the symmetric warm two stream instability, i.e., the electron distribution function in the VP system starts with the following unstable initial condition, 2 v2 f (x, v, t = 0) = √ (1 + 5v 2 )(1 + α((cos(2kx) + cos(3kx))/1.2 + cos(kx)) exp(− ), (6.9) 2 7 2π where α = 0.01, k = 0.5 and the length of the domain in the x direction is L =

2π k

and

in the v direction vmax = 5. The background ion distribution function is fixed, uniform and chosen so that the total net charge density for the system is zero. Figure 6.3 shows 27

numerical solutions of phase space profiles at T = 53 from the MPP SL DG scheme with different solution spaces Vhk . Consistent numerical results are generated for all Vhk compared with those presented in [14, 25]. The higher the k is, the better resolution of the solution is shown. Numerical solution from MPP SL DG scheme with P 1 , but with much more refined mesh is also shown as the reference solution in Figure 6.3, indicating the convergence of the numerical solution. To see more details, we also provide the contour plots with 30 equally spaced contour curves on [0, 1] in Figure 6.4. Numerical results are observed to be consistent with each other with some small but artificial numerical oscillations (see the bottom plots in Figure 6.4). Figure 6.5 shows the time development of the discrete L1 norm, L2 norm, kinetic energy and entropy for MPP SL DG scheme with different solution spaces Vhk . It is observed that all the schemes are conservative in the L1 norm, because of the conservative and MPP nature of the numerical schemes. Schemes with higher degree in polynomial solution space in general do a better job in preserving other relevant physical norms. Example 6.4. (Two stream instability [11]) Consider the symmetric two stream instability, similar as in [25], f (x, v, t = 0) =

1 √

2vth

     (v − u)2 v+u + exp − 2 (1 + 0.05 cos(kx)) exp − 2 2vth 2vth 2π

with u = 0.99, vth = 0.3 and k =

2 . 13

(6.10)

The background ion distribution function is fixed,

uniform and chosen so that the total net charge density for the system is zero. Our numerical simulation parameters are vmax = 5, Nx = 512, Nv = 512, CF L = 3 for all schemes. Figure 6.6 shows numerical solutions of phase space profiles at T = 70 from MPP SL DG schemes with different solution spaces Vhk and with two sets of different meshes. To see more details, we also provide the contour plots with 30 equally spaced contour curves on [0, 1] in Figure 6.7. Figure 6.8 shows the time development of the discrete L1 norm, L2 norm, kinetic energy and entropy for MPP SL DG scheme with different Vhk . Again all schemes are conservative in the L1 norm, because of the conservative and MPP nature of the numerical schemes. Schemes with higher degree in polynomial solution space in general do a better job in preserving other relevant physical norms.

28

7

Conclusions

In this paper, we propose a maximum principle preserving (MPP) or positivity preserving (PP) semi-Lagrangian (SL) discontinuous Galerkin (DG) method for scalar linear advection equation with general variable coefficients. The scheme inherits the flexibility, compactness and hp-adaptivity capability from the DG formulation, as well as the allowance of large time step evolution in the SL framework. A MPP/PP limiter is applied to ensure MPP/PP property, when the exact solution of the PDE has such properties, without affecting the original high order accuracy. We analyze the stability and accuracy properties of the scheme by establishing the connection with the weak form of the Lagrangian Galerkin method [22]. We apply the method to several test examples, such as advection in incompressible flows and the VP system. Consistent numerical results are observed with those in the literature.

References [1] D. Arnold, F. Brezzi, B. Cockburn, and L. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis, 39 (2002), pp. 1749–1779. [2] M. Begue, A. Ghizzo, P. Bertrand, E. Sonnendrucker, and O. Coulaud, Two-dimensional semi-Lagrangian Vlasov simulations of laser–plasma interaction in the relativistic regime, Journal of Plasma Physics, 62 (1999), pp. 367–388. [3] N. Besse and E. Sonnendrucker, Semi-Lagrangian schemes for the Vlasov equation on an unstructured mesh of phase space, Journal of Computational Physics, 191 (2003), pp. 341–376. [4] J. A. Carrillo and F. Vecil, Nonoscillatory interpolation methods applied to Vlasov-based models, SIAM Journal on Scientific Computing, 29 (2007), pp. 1179–1206. [5] C. Cheng and G. Knorr, The integration of the Vlasov equation in configuration space, Journal of Computational Physics, 22 (1976), pp. 330–351.

29

[6] B. Cockburn, C. Johnson, C.-W. Shu, and E. Tadmor, Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Springer, New York, 1998. [7] B. Cockburn, S. Lin, and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one-dimensional systems, Journal of Computational Physics, 84 (1989), pp. 90–113. [8] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, (1989), pp. 411–435. , The Runge-Kutta discontinuous Galerkin method for conservation laws V: multi-

[9]

dimensional systems, Journal of Computational Physics, 141 (1998), pp. 199–224. [10]

, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, Journal of Scientific Computing, 16 (2001), pp. 173–261.

[11] N. Crouseilles, M. Mehrenberger, and E. Sonnendrucker, Conservative semi-Lagrangian schemes for Vlasov equations, Journal of Computational Physics, 229 (2010), pp. 1927–1953. [12] J. Douglas Jr and T. Russell, Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures, SIAM Journal on Numerical Analysis, 19 (1982), pp. 871– 885. [13] F. Filbet and E. Sonnendrucker, Comparison of eulerian Vlasov solvers, Computer Physics Communications, 150 (2003), pp. 247–266. ¨cker, Comparison of Eulerian Vlasov solvers, Com[14] F. Filbet and E. Sonnendru puter Physics Communications, 150 (2003), pp. 247–266. ¨cker, and P. Bertrand, Conservative numerical [15] F. Filbet, E. Sonnendru schemes for the Vlasov equation, Journal of Computational Physics, 172 (2001), pp. 166– 187. 30

[16] S. Gottlieb, D. Ketcheson, and C.-W. Shu, High order strong stability preserving time discretizations, Journal of Scientific Computing, 38 (2009), pp. 251–289. [17] R. E. Heath, I. Gamba, P. Morrison, and C. Michler, A discontinuous Galerkin method for the Vlasov-Poisson system, Numerische Mathematik, 53 (1988), pp. 459–483. [18] F. Huot, A. Ghizzo, P. Bertrand, E. Sonnendrucker, and O. Coulaud, Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov– Maxwell system, Journal of Computational Physics, 185 (2003), pp. 512–531. [19] G.-S. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes, Journal of Computational Physics, 126 (1996), pp. 202–228. [20] E. Kubatko, C. Dawson, and J. Westerink, Time step restrictions for RungeKutta discontinuous Galerkin methods on triangular grids, Journal of Computational Physics, 227 (2008), pp. 9697–9710. [21] R. LeVeque, High-resolution conservative algorithms for advection in incompressible flow, SIAM Journal on Numerical Analysis, (1996), pp. 627–665. [22] K. Morton, A. Priestley, and E. Suli, Stability of the Lagrange-Galerkin method with non-exact integration, Mod´elisation Math´ematique et Analyse Num´erique, 22 (1988), pp. 625–653. [23] J.-M. Qiu and A. Christlieb, A Conservative high order semi-Lagrangian WENO method for the Vlasov Equation, Journal of Computational Physics, 229 (2010), pp. 1130–1149. [24] J.-M. Qiu and C.-W. Shu, Conservative high order semi-Lagrangian finite difference WENO methods for advection in incompressible flow, Journal of Computational Physics, accepted. [25]

, Conservative semi-Lagrangian finite difference WENO formulations with applications to the Vlasov equation, Communications in Computational Physics, (accepted).

31

[26] M. Restelli, L. Bonaventura, and R. Sacco, A semi-Lagrangian discontinuous Galerkin method for scalar advection by incompressible flows, Journal of Computational Physics, 216 (2006), pp. 195–215. [27] J. Rossmanith and D. Seal, A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations, Journal of Computational Physics, submitted. [28] E. Sonnendruecker, J. Roche, P. Bertrand, and A. Ghizzo, The semiLagrangian method for the numerical resolution of the Vlasov equation, Journal of Computational Physics, 149 (1999), pp. 201–220. [29] A. Staniforth and J. Cote, Semi-Lagrangian integration schemes for atmospheric modelsA review, Monthly Weather Review, 119 (1991), pp. 2206–2223. [30] E. Suli, Convergence and nonlinear stability of the Lagrange-Galerkin method for the Navier-Stokes equations, Numerische Mathematik, 53 (1988), pp. 459–483. [31] T. Umeda, M. Ashour-Abdalla, and D. Schriver, Comparison of numerical interpolation schemes for one-dimensional electrostatic Vlasov code, Journal of Plasma Physics, 72 (2006), pp. 1057–1060. [32] H. Wang, R. Ewing, G. Qin, S. Lyons, M. Al-Lawatia, and S. Man, A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations, Journal of Computational Physics, 152 (1999), pp. 120–163. [33] T. Yabe, F. Xiao, and T. Utsumi, The constrained interpolation profile method for multiphase analysis, Journal of Computational Physics, 169 (2001), pp. 556–593. [34] X. Zhang and C.-W. Shu, A genuinely high order total variation diminishing scheme for one-dimensional scalar conservation laws, SIAM Journal on Numerical Analysis, 48 (2010), pp. 772–795. [35] X. Zhang and C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, Journal of Computational Physics, 229 (2010), pp. 3091–3120. 32

[36]

, On positivity preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, Journal of Computational Physics, 229 (2010), pp. 8918–8934.

[37] T. Zhou, Y. Guo, and C.-W. Shu, Numerical study on Landau damping, Physica D: Nonlinear Phenomena, 157 (2001), pp. 322–333.

33

Figure 5.4: Plots of the numerical solution of MPP SL DG with different solution space Vhk for equation (5.4) with CF L = 2.2 at T = 2π. The numerical mesh is uniform and is 192 × 192 for P 0 solution space, 96 × 96 for P 1 solution space, 64 × 64 for P 2 solution space and 48 × 48 for P 3 solution space. Different numerical meshes are chosen for different Vhk , so that the overall degrees of freedom are the same for different solution spaces.

34

Figure 5.5: Plots of slides of numerical solutions for equation (5.4) at Y = π/2 (top left), X = π (top right) and Y = 3π/2 (bottom left) in Figure 5.4. The bottom right figure is the zoom-in plot for the bottom left figure around the mid upper region.

35

Figure 5.6: Plots of the numerical solution of Strang split PP SL DG for equation (5.5) with different solution space Vhk , CF L = 0.2 at T = 1.5. The numerical mesh is uniform and is 96 × 96 for different solution spaces.

36

Figure 5.7: Plots of slides of numerical solution in Figure 5.6 for equation (5.4) at Y = π/2 (top), X = π (middle) and Y = 3π/2 (bottom) with CF L = 0.2 at T = 1.5.

37

Figure 5.8: Plots of the numerical solution of PP SL DG scheme with P 1 solution space for equation (5.5) with CF L = 0.2, T = 1.5 and final integration time 0.75, when the initial profile is quite deformed. The numerical mesh is 96 × 96.

38

Figure 5.9: Plots of the numerical solution of Strang split MPP SL DG for equation (5.5) with different solution space Vhk , CF L = 0.2 at T = 1.5. The numerical mesh is uniform and is 96 × 96 for different solution spaces.

39

Figure 5.10: Plots of slides of numerical solution in Figure 5.9 for equation (5.4) at Y = π/2 (top), X = π (middle) and Y = 3π/2 (bottom) with CF L = 0.2 at T = 1.5.

40

Figure 5.11: Plots of the numerical solution of MPP SL DG scheme with P 1 solution space for equation (5.5) with CF L = 0.2, T = 1.5 and final integration time 0.75, when the initial profile is quite deformed. The numerical mesh is 96 × 96.

41

Figure 6.1: Weak Landau damping: time evolution of electric field in L2 norm for P 1 (top left), P 2 (top middle) and P 3 (top right), L1 (mid left) and L2 (mid right) norms of the solution as well as the discrete kinetic energy (lower left) and entropy (lower right).

42

Figure 6.2: Strong Landau damping: time evolution of electric field in L2 (upper) norm, L1 (mid left) and L2 (mid right) norms of the solution as well as the discrete kinetic energy (lower left) and entropy (lower right). 43

Figure 6.3: Phase space plots of the two stream instability at T = 53 using MPP SL DG with P 1 (upper left), P 2 (upper right), P 3 (lower left) with mesh 64 × 128 and CF L = 3 for all test cases. The lower right figure is a reference solution produced by the scheme with P 1 solution space but with refined mesh 256 × 512.

44

Figure 6.4: Phase space plots of the two stream instability at T = 53 using MPP SL DG with P 1 (upper left), P 2 (upper right), P 3 (mid left) with mesh 64 × 128 and CF L = 3 for all test cases. The mid right figure is a reference solution produced by the scheme with P 1 solution space but with refined mesh 256 × 512. The two plots in the bottom row are the corresponding zoom in plots for the two plots in the middle row. 45

Figure 6.5: Two-stream instability: time evolution of the L1 (upper left) and L2 (upper right) norms of the solution as well as the discrete kinetic energy (lower left) and entropy (lower right).

46

Figure 6.6: Phase space plots of the two stream instability at T = 70 using MPP SL DG with P 1 (upper row), P 2 (middle row), P 3 (lower row) with mesh 64 × 128 (left) and 128 × 256 (right). CF L = 3 for all test cases.

47

Figure 6.7: Phase space plots of the two stream instability at T = 70 using MPP SL DG with P 1 (upper left), P 2 (upper right), P 3 (lower left) with mesh 128 × 256, CF L = 3 for all test cases. The lower right plot is the zoom in plot of the lower left plot.

48

Figure 6.8: Two-stream instability: time evolution of the L1 (upper left) and L2 (upper right) norms of the solution as well as the discrete kinetic energy (lower left) and entropy (lower right).

49