An Eulerian-Lagrangian discontinuous Galerkin method for transport problems and its application to nonlinear dynamics
AAN EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKINMETHOD FOR TRANSPORT PROBLEMS ANDITS APPLICATION TO NONLINEAR DYNAMICS ∗ XIAOFENG CAI † , JING-MEI QIU ‡ , AND
YANG YANG § Abstract.
We propose a new Eulerian-Lagrangian (EL) discontinuous Galerkin (DG) method.The method is designed as a generalization of the semi-Lagrangian (SL) DG method for linearadvection problems proposed in [J. Sci. Comput. 73: 514-542, 2017], which is formulated based onan adjoint problem and tracing upstream cells by tracking characteristics curves highly accurately.In the SLDG method, depending on the velocity field, upstream cells could be of arbitrary shape.Thus, a more sophisticated approximation to sides of the upstream cells is required to get high orderapproximation. For example, quadratic-curved (QC) quadrilaterals were proposed to approximateupstream cells for a third-order spatial accuracy in a swirling deformation example. In this paper, forlinear advection problems, we propose a more general formulation, named the ELDG method. Thescheme is formulated based on a modified adjoint problem for which the upstream cells are alwaysquadrilaterals, which avoids the need to use QC quadrilaterals in the SLDG algorithm. The newlyproposed ELDG method can be viewed as a new general framework, in which both the classicalEulerian Runge-Kutta DG formulation and the SL DG formulation can fit in. Numerical results onlinear transport problems, as well as the nonlinear Vlasov and incompressible Euler dynamics usingthe exponential RK time integrators, are presented to demonstrate the effectiveness of the ELDGmethod.
Key words.
Eulerian-Lagrangian; discontinuous Galerkin; mass conservative; semi-Lagrangian;Vlasov simulations; characteristics.
AMS subject classifications.
1. Introduction.
We propose a new Eulerian-Lagrangian (EL) discontinuousGalerkin (DG) method for a model transport equation in the form of(1.1) u t + ∇ · ( P ( u ; x , t ) u ) = 0 , ( x , t ) ∈ R d × [0 , T ] , which could come from a wide range of application fields including fluid dynamics,climate modeling, and kinetic description of plasma. There are three main classes ofcomputational methods for solving (1.1): Lagrangian, Eulerian and semi-Lagrangian(SL). Each class of methods has their own advantages and limitations. The Lagrangianmethod is particle based, works efficiently for high dimensional problems, but suffersfrom statistical noises; while the latter two methods are mesh-based method, can bedesigned to be of high order accurate, but suffers from the curse of dimensionality.The main difference between Eulerian and SL methods is the space-time region in con-sideration: the Eulerian method performs numerical discretizations with fixed spatiallocations in time; while the semi-Lagrangian method usually do that along convection ∗ Funding:
Research of the first and second author is supported by NSF grant NSF-DMS-1818924,Air Force Office of Scientific Computing FA9550-18-1-0257 and University of Delaware. Research ofthe third author is supported by NSF grant DMS-1818467. † Department of Mathematical Sciences, University of Delaware, Newark, DE, 19716, USA. ([email protected]). ‡ Department of Mathematical Sciences, University of Delaware, Newark, DE, 19716, USA.([email protected]). § Department of Mathematical Sciences, Michigan Technological University, Houghton, MI 49931,USA. ([email protected]). 1 a r X i v : . [ m a t h . NA ] F e b X. CAI, J.-M. QIU, AND Y. YANG characteristics. When characteristics are tracked accurately, semi-Lagrangian meth-ods often allow much larger time stepping sizes than their Eulerian counterparts.Among different classes of SL methods in the literature, we would like to mention afew closely related ones that are developed in the finite element framework. There is aline of research work along Eulerian Lagrangian Localized Adjoint Methods (ELLAM)[8]. ELLAM introduces an adjoint problem for the test function in the continuousfinite element framework and has a broad range of influence in different applicationfronts [32, 29]. Compared with ELLAM, the SLDG [5] is being developed in thediscontinuous Galerkin finite element framework. SL schemes could be developedbase on forward [3] or backward characteristics tracing. Here we choose to developour schemes base on backward characteristics tracing.In this paper, we propose a new ELDG method that is mesh-based, and is a gen-eralized framework of the SL DG method developed earlier [5]. It is designed to takeadvantage of information propagation along characteristics as in a SL method, andmaintain essential properties of the SLDG method on mass conservation, high orderspatial and temporal accuracy, and allowing for extra large time steps with stability.We first focus on developing the ELDG algorithm for linear transport problems. Anew ingredient of the method is the introduction of a modified adjoint problem forthe test function.
The velocity field of the modified adjoint problem is a linear func-tion that approximates that of the original transport problem.
There are two positiveconsequences of such modification. One is that the test function remains in the same P k polynomial spaces, whereas in the SLDG setting the test function does not neces-sarily remain in P k and needs to be approximated. In fact, a close connection can bedrawn between the ELDG method and the Arbitrary Lagrangian Eulerian (ALE) DGmethod [23], when we view the space-time region in the ELDG method as a dynamicmoving mesh. The second advantage brought by the modified adjoint problem is thatthe shape of upstream cells is always quadrilaterals in a 2D setting. For a generalvariable coefficient problem, upstream cells of the SLDG method could be of arbitraryshape and needs to be better approximated. In [5], we propose to use quadratic curvesin approximating sides of upstream cells, so that we have third order spatial accuracy.Such a practice is difficult be further generalized to schemes with even higher orderaccuracy, and for problems in higher-dimensions. With the newly ELDG method, nocurves are needed to better approximate upstream cells. A direct generalization ofthe algorithm to higher dimensional problems can be similarly done in principle.Due to the approximate nature of the velocity field in the modified adjoint prob-lem, there is an extra flux term taking account of the difference between velocityfields from the modified adjoint problem and the original problem. The newly pro-posed ELDG scheme evolves this extra flux term in a similar spirit to the classicalEulerian RKDG method [12]. The ELDG scheme is designed base on the integralform of the equation over characteristics-related space-time regions; yet we transformsuch integral formulation into a time-differential form, for which the method-of-linesstrong-stability preserving (SSP) Runge-Kutta (RK) can be directly applied. Here,we would like to mention the Eulerian Lagrangian weighted essentially non-oscillatoryschemes developed in [20, 21], for which a different way of treating time integrationis proposed.As nonlinear applications of the ELDG algorithm, we consider the nonlinearVlasov-Poisson system, the guiding center Vlasov model as well as the incompress-ible Euler equations. Here, we couple the ELDG algorithm with the RK exponentialintegrator [9, 4] to realize a uniformly high order spatial-temporal discretization ofnonlinear transport. In particular, the RK exponential integrator decomposes a time N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD
2. ELDG formulation for 1D linear transport problems.
To illustrate thekey idea of the ELDG scheme, we start from a 1D linear transport equation in thefollowing form(2.1) u t + ( a ( x, t ) u ) x = 0 , x ∈ [ x a , x b ] . For simplicity, we assume periodic boundary conditions, and the velocity field a ( x, t ) isa continuous function of space and time. We perform a partition of the computationaldomain x a = x < x < · · · < x N + = x b . Let I j = [ x j − , x j + ] denote an elementof length ∆ x j = x j + − x j − and define ∆ x = max j ∆ x j . We define the finitedimensional approximation space, V kh = { v h : v h | I j ∈ P k ( I j ) } , where P k ( I j ) denotesthe set of polynomials of degree at most k . For this finite-dimensional space, weintroduce a set of basis functions { Ψ j,m ( x ) } ≤ j ≤ N, ≤ m ≤ k . We also introduce a set ofbasis functions { ψ j,m ( x, t ) } ≤ j ≤ N, ≤ m ≤ k , which will be used in an adjoint problem.The subscripts of Ψ j,m ( x ) and ψ j,m ( x, t ) are often omitted, when there is no risk ofambiguity. Moreover, we define t n to be the n − th time level, and ∆ t = t n +1 − t n tobe the time-stepping size. The SLDG method proposed in [5] isformulated based on an adjoint problem of (2.1) with ∀ Ψ ∈ P k ( I j ),(2.2) (cid:40) ψ t + a ( x, t ) ψ x = 0 , t ∈ [ t n , t n +1 ] ,ψ ( t = t n +1 ) = Ψ( x ) , for which the solution ψ stays constant along characteristic trajectories. It was shownin [19] that(2.3) ddt (cid:90) (cid:101) I j ( t ) u ( x, t ) ψ ( x, t ) dx = 0 , where (cid:101) I j ( t ) is a dynamic interval bounded by characteristics emanating from cellboundaries of I j at t = t n +1 , see Figure 1 for illustration. An SL time discretizationof (2.3) leads to(2.4) (cid:90) I j u n +1 Ψ dx = (cid:90) I (cid:63)j u ( x, t n ) ψ ( x, t n ) dx, where I (cid:63)j = [ x (cid:63)j − , x (cid:63)j + ] with x (cid:63)j ± = (cid:101) x j ± ( t n ) being the foots of trajectory at t n emanating from ( x j ± , t n +1 ). In order to update the numerical solution u n +1 , we varythe test function Ψ as basis of V kh and evaluate the right-hand side (RHS) integral of(2.4) properly. The detailed procedures can be found in [5]. X. CAI, J.-M. QIU, AND Y. YANG t n t n +1 x j − x j + x (cid:63)j − x (cid:63)j + (cid:101) I j ( t n ) = I (cid:63)j (cid:101) I j ( t n +1 ) = I j (cid:101) I j ( t ) d ˜ x ( t ) dt = a (˜ x ( t ) , t )where ˜ x ( tn +1) = xj − Fig. 1 . Illustration for the space-time region for the SLDG formulation.
The newly proposed ELDG method differsfrom the SLDG method [5] in the formulation of a modified adjoint problem for the testfunction ψ . To introduce the scheme, we first introduce the modified adjoint problemand the associated space-time region Ω j ; then we derive a semi-discrete version ofthe ELDG scheme based on the space-time region of Ω j ; finally a method-of-timesRunge-Kutta method is applied for time marching. (1) A modified adjoint problem. We consider the adjoint problem with ∀ Ψ ∈ P k ( I j ) on the time interval [ t n , t n +1 ]:(2.5) (cid:40) ψ t + α ( x, t ) ψ x = 0 , t ∈ [ t n , t n +1 ] ,ψ ( t = t n +1 ) = Ψ( x ) , with α ( x, t ) being a bilinear function of ( x, t ) designed by three steps below:1. On I j at t n +1 : we let α ( x, t n +1 ) be a linear polynomial on I j interpolating a ( x, t n +1 ) at cell boundaries,(2.6) α ( x j ± , t n +1 ) = a ( x j ± , t n +1 ) . = ν j ± . That is,(2.7) α ( x, t n +1 ) = − ν j − x − x j + ∆ x j + ν j + x − x j − ∆ x j ∈ P ( I j ) .
2. We define a space-time region Ω j = ˜ I j ( t ) × [ t n , t n +1 ] with the dynamic in-terval, ˜ I j ( t ) = [˜ x j − ( t ) , ˜ x j + ( t )], t ∈ [ t n , t n +1 ], where ˜ x j ± ( t ) = x j ± +( t − t n +1 ) a ( x j ± , t n +1 ) emanating from cell boundaries x j ± with slopes a ( x j ± , t n +1 ). It will become clear after the third step that the space-timeregion Ω j is the dynamic characteristic region of the modified adjoint problem(2.5). We let I (cid:63)j . = ˜ I j ( t n ) be the upstream cell of I j at t n . See the left panelin Figure 2 for illustration.3. On ˜ I j ( t ) for [ t n , t n +1 ) : let ˜ x ( t ; ( ξ, t n +1 )) be a straight line emanating fromany point ξ ∈ I j at t n +1 and with the slope α ( ξ, t n +1 ). That is,(2.8) ddt ˜ x ( t ; ( ξ, t n +1 )) = α ( ξ, t n +1 ) , ˜ x ( t n +1 ; ( ξ, t n +1 )) = ξ. Then(2.9) ˜ x ( τ ; ( ξ, t n +1 )) = ξ − α ( ξ, t n +1 )( t n +1 − τ ) , ∀ τ ∈ [ t n , t n +1 ) . N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD α (˜ x ( τ ; ( ξ, t n +1 )) , τ ) = α ( ξ, t n +1 ) , τ ∈ [ t n , t n +1 ) . We would like to point out a few facts about Ω j and the modified adjoint problem(2.5): • From the construction of Ω j and α ( x, t ) of the modified adjoint problem(2.5), it can be easily checked that, (2.8) is the characteristics equation forthe modified adjoint problem (2.5). • ˜ x ( τ ; ( ξ, t n +1 )) satisfying eq. (2.9) is a linear function of ξ and τ ; the Jacobianis(2.11) ∂ ˜ x ( τ ; ( ξ, t n +1 )) ∂ξ = 1 − ν j + − ν j − ∆ x j ( t n +1 − τ ) , which will become useful later in implementation. In particular, ∂ ˜ x ( t n ; ( ξ, t n +1 )) ∂ξ = 1 − ∆ t ν j + − ν j − ∆ x j . • In order for the characteristics not crossing each other, one has to enforce thecondition of ∂ ˜ x ( t n ;( ξ,t n +1 )) ∂ξ ≥
0, which implies the time step constraint(2.12) ∆ t ≤ min j ∆ x j max( ν j + − ν j − , . • For the modified adjoint problem, the solution ψ stays constant along char-acteristics (2.5), therefore we have(2.13) ψ (˜ x ( τ ; ( ξ, t n +1 )) , τ ) = Ψ( ξ ) ∈ P k ( I j ) , ∀ τ ∈ [ t n , t n +1 ] . If we consider a transformation between x ∈ ˜ I j to a reference interval ξ ∈ I j ,see Figure 2, eq. (2.13) indicates that the test function ψ (˜ x ( τ ; ( ξ, t n +1 )) , τ )in the ξ coordinate remains the same as the classical test function Ψ( ξ ), i.e.standard basis functions in P k ( I j ). Ω j t n t n +1 x j − x j + x (cid:63)j − x (cid:63)j + (cid:101) I j ( t n ) = I (cid:63)j (cid:101) I j ( t n +1 ) = I j (cid:101) I j ( t ) ν j − ν j + ξ = x j − / ξ = x j +1 / Fig. 2 . Illustration for the mapping between dynamic element (cid:101) I j ( t ) (left) and the iso-parametricelement (right). X. CAI, J.-M. QIU, AND Y. YANG (2) Formulation of the semi-discrete ELDG scheme.
In order to formulate thescheme, we integrate (2.1) · ψ + (2.5) · u over Ω j , which gives the following identity,(2.14) (cid:90) Ω j [(2.1) · ψ + (2.5) · u ] dxdt = 0 . That is,0 = (cid:90) t n +1 t n (cid:90) ˜ I j ( t ) ( u t ψ + uψ t ) dxdt + (cid:90) t n +1 t n (cid:90) ˜ I j ( t ) (( a ( x, t ) u ) x ψ + α ( x, t ) ψ x u ) dxdt = (cid:90) t n +1 t n (cid:90) ˜ I j ( t ) ( uψ ) t dxdt + (cid:90) t n +1 t n (cid:90) ˜ I j ( t ) (( auψ ) x − auψ x + αψ x u ) dxdt = (cid:90) t n +1 t n (cid:34) ddt (cid:90) ˜ I j ( t ) uψdx − αuψ (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t )˜ x j − ( t ) + auψ (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t )˜ x j − ( t ) + (cid:90) ˜ I j ( t ) ( α − a ) uψ x dx (cid:35) dt = (cid:90) t n +1 t n (cid:34) ddt (cid:90) ˜ I j ( t ) uψdx + ( a − α ) uψ (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t )˜ x j − ( t ) − (cid:90) ˜ I j ( t ) ( a − α ) uψ x dx (cid:35) dt. (2.15)Letting F ( u ) . = ( a − α ) u , the time differential form of (2.15) gives(2.16) ddt (cid:90) ˜ I j ( t ) ( uψ ) dx = − ( F ψ ) (cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + ( F ψ ) (cid:12)(cid:12)(cid:12) ˜ x j − ( t ) + (cid:90) ˜ I j ( t ) F ψ x dx. Notice that the dynamic interval of ˜ I j ( t ) can always be linearly mapped to a referencecell ξ ∈ I j , see the right plot in Figure 2, then eq. (2.16) in the ξ -coordinate becomes(2.17) ddt (cid:90) I j ( u Ψ( ξ )) ∂ ˜ x ( t ; ( ξ, t n +1 )) ∂ξ dξ = − ( F Ψ) (cid:12)(cid:12)(cid:12) ξ = x j + 12 + ( F Ψ) (cid:12)(cid:12)(cid:12) ξ = x j − + (cid:90) I j F Ψ ξ dξ. The DG discretization [13, 12] of (2.17) is to find u h ( ξ, t ) ∈ P k ( I j ) as the approximatesolution of u (˜ x ( t ; ( ξ, t n +1 )) , t ) on ˜ I j ( t ), so that for ∀ Ψ ∈ P k ( I j ),(2.18) ddt (cid:90) I j u h Ψ ∂ ˜ x ( t ; ( ξ, t n +1 )) ∂ξ dξ = − ˆ F j + Ψ( x − j + ) + ˆ F j − Ψ( x + j − ) + (cid:90) I j F Ψ ξ dξ. Notice here u h could be discontinuous across x (cid:63)j − . In this paper, we choose ˆ F as amonotone flux, e.g. the Lax-Friedrichs flux(2.19) ˆ F ( u − , u + ) = 12 ( F ( u − ) + F ( u + )) − α u + − u − ) , α = max u | F (cid:48) ( u ) | ;and we use Gauss quadrature rules with k + 1 quadrature points to approximate theintegral term (cid:82) I j F ( u h )Ψ ξ dξ on the RHS of the equation (2.18). (3) RK time discretization and fully discrete scheme. We can write the semi-discrete scheme (2.18) into a form of ordinary differential equations (ODEs) with aninitial condition. We let ˜ U ( t ) be a vector in R N ( k +1) which consists of degrees offreedom { (cid:82) ˜ I j ( t ) u h ( x, t ) ψ j,m ( x, t ) dx . = ˜ U j,m ( t ) } ≤ j ≤ N, ≤ m ≤ k , and denote the spatial N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD L (cid:16) ˜ U ( t ) , t (cid:17) . Then the semi-discretescheme (2.18) can be written as ∂∂t ˜ U ( t ) = L (cid:16) ˜ U ( t ) , t (cid:17) , ˜ U ( t n ) = ˜ U n . (2.20)There are two main steps involved here.1. Obtain the initial condition of (2.18) by an L projection of u h on up-stream cells ˜ I j by SLDG method. In particular, ˜ U n consists of the numericalsolutions ˜ U nj,m of the SLDG scheme [5] for approximating (cid:90) ˜ I j ( t n ) u h ( x, t n ) ψ j,m ( x, t n ) dx. Update (2.20) from ˜ U n to ˜ U n +1 . we apply the SSP explicit RK meth-ods [31] as in a method-of-lines approach. In particular, the time-marchingalgorithm using an s -stage RK method follows the procedure below:(a) Get the mesh information of the dynamic element ˜ I ( l ) j , l = 0 , · · · , s onRK stages by eq. (2.9).(b) For RK stages i = 1 , · · · , s , compute˜ U ( i ) = i − (cid:88) l =0 (cid:104) α il ˜ U ( l ) + β il ∆ t n L (cid:16) ˜ U ( l ) , t n + d l ∆ t n (cid:17)(cid:105) , (2.21)where α il and β il are related to RK methods. They are provided inTable 1 for the second order and third order SSP RK methods.Note that ˜ U n is evaluated by the SLDG scheme in x -coordinate, while ˜ U ( i ) in the each time stage is updated with respect to the reference ξ coordinate. Table 1
Parameters of some practical Runge-Kutta time discretizations.
Order α il β il d l
12 12
13 1 1 0
34 14
23 12
Theorem (Mass conservation) Given a DG solution u h ( x, t n ) ∈ V kh andassuming the boundary condition is periodic, the proposed fully discrete ELDG schemewith SSP RK time discretization of (2.20) is locally mass conservative. In particular, N (cid:88) i =1 (cid:90) I j u h ( x, t n +1 ) dx = N (cid:88) i =1 (cid:90) I j u h ( x, t n ) dx. Proof.
It can be proved by letting ψ = 1, the conservative form of integrating F function with unique flux at cell boundaries, as the mass conservation property ofSLDG scheme [5]. We skip details for brevity. X. CAI, J.-M. QIU, AND Y. YANG
A few remarks are in order for the proposed ELDG scheme, in comparison withexisting SLDG [5], RKDG [12] and ALE DG [23] methods in the literature. Theseremarks also apply to the 2D ELDG scheme in the next section.
Remark α ( x, t ) is an approximation of a ( x, t ). While the characteristics induced by a ( x, t ) could be curves and the test function φ satisfying eq. (2.2) may no longerbe polynomials, the characteristics induced by α ( x, t ) are straight lines and the testfunction φ remains a P k polynomial on ˜ I j ( t ). The difference, between α ( x, t ) andexact slopes a ( x, t ) for characteristic curves, is taken into account by the F functionin (2.18). Remark α = a , the ELDG method becomes the SLDG method [5] andthe scheme is unconditionally stable. In the special case of α ( x, t ) = 0, the ELDGmethod becomes the classical RKDG method [15]. In the general setting that α ap-proximates (but not exactly equals) a , the ELDG method enables larger time stepconstraint for stability than the classical DG scheme. One can compare the time stepconstraint (2.23) to that of a classical Eulerian DG method. Remark I j at t n +1 and the upstream cells I (cid:63)j at t n in a movingmesh setting, the formulation of ELDG (2.18) is the same as the ALE DG method[23] and the quasi-Lagrangian moving mesh discontinuous Galerkin method [25]. Afundamental difference between the ELDG and ALE DG methods is that the latterone is formulated based on a set of moving mesh, whereas the ELDG method inthis paper is based on a fixed set of mesh. As a result, the ELDG method avoid thecomplication of mesh distortion as in an ALE DG method. In fact, the ELDG methodcan be viewed as a combination of SLDG algorithm in evaluating ˜ U n and an ALEDG method in updating solutions from ˜ U n to ˜ U n +1 . Remark F = ( a − α ) u , thus an empirical time step stabilityconstraint of the proposed ELDG method is(2.22) ∆ t ≤ ∆ x (2 k + 1) max | a ( x, t ) − α ( x, t ) | , with k being the polynomial degree of the DG method. Combine this with (2.12)gives(2.23) ∆ t ≤ ∆ x max { (2 k + 1) max | a ( x, t ) − α ( x, t ) | , a ( x j + , t n +1 ) − a ( x j − , t n +1 ) } . For a smooth function a , from the construction of α function as previously describedand by Taylor expansions, we have α − a = O (∆ t )+ O (∆ x ). Combining this estimatewith (2.23) give the time step constraint for stability of ELDG∆ t ∼ ∆ x . This is consistent with our numerical observations presented in Section 5.
N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD Remark u t + u x =0 with α ( x, t ) for the adjoint problem being a constant α (cid:54) = 1 could be obtained by thestability of an L projection as in an SLDG scheme [27], together with the stabilityof a fully discrete ALE DG method [37]. A rigorous analysis is subject to furtherinvestigation. Remark α ( x, t n +1 ) is constructed asa linear function interpolating a ( x, t ) at cell boundaries. Alternatively, for (2.1),one can track characteristics from cell boundaries at t n +1 , i.e. from ( x j ± / , t n +1 )find their characteristics feet ( x (cid:63)j ± / , t n +1 ). Then α ( x j ± / , t n +1 ) can be obtainedas the slope of the straight time connecting ( x j ± / , t n +1 ) and ( x (cid:63)j ± / , t n +1 ), i.e. α ( x j ± / , t n +1 ) = x j ± / − x (cid:63)j ± / ∆ t . We name the ELDG scheme with such constructionof α function as ‘ELDG-ST2’, and the ELDG scheme with α ( x, t ) defined by eq. (2.7)and (2.10) as ‘ELDG-ST1’ in later parts of this paper.
3. The ELDG algorithm for 2D transport problems..
The design of the2D ELDG algorithm shares a similar spirit as the 1D case. We consider a lineartransport equation(3.1) u t + ( a ( x, y, t ) u ) x + ( b ( x, y, t ) u ) y = 0 , ( x, y ) ∈ Ω . For simplicity, we assume the computational domain Ω is rectangular, boundary con-ditions are periodic, and the velocity field ( a ( x, y, t ) , b ( x, y, t )) is a continuous functionof space and time. We partition the domain Ω by a set of non-overlapping rectangularelements A j , j = 1 , · · · , J , and define the finite dimensional DG approximation space, V kh = { v h : v h | A j ∈ P k ( A j ) } , where P k ( A j ) denotes the set of polynomials of degreeat most k over A j = [ x lj , x rj ] × [ y bj , y tj ] with element center (cid:18) x j = x lj + x rj , y j = y bj + y tj (cid:19) and sizes, ∆ x j = x rj − x lj , ∆ y j = y tj − y bj . Let n k be the dimension of P k ( A j ). (1) A modified adjoint problem for the 2D transport problem. To derive a2D ELDG formulation, we consider a modified adjoint problem at ˜ A j ( t ) on the timeinterval t ∈ [ t n , t n +1 ]:(3.2) ψ t + α ( x, y, t ) ψ x + β ( x, y, t ) ψ y = 0 , ψ ( x, y, t = t n +1 ) = Ψ( x, y ) ∈ P k ( A j ) , where ( α, β ) are bilinear functions on A j at t n +1 defined as described below. Notation-wise, we let ˜ A j ( t ) , t ∈ [ t n , t n +1 ] be the dynamic characteristic element of the modifiedadjoint problem (3.2) with (˜ x ( t ) , ˜ y ( t )) ∈ ˜ A j ( t ) that satisfies (3.4) emanating from( x, y ) of A j at t n +1 . We also let A (cid:63)j . = ˜ A j ( t n ) be the upstream cell of A j at t n and letΩ j be the region of which ( x, y, t ) ∈ ˜ A j ( t ) × [ t n , t n +1 ].1. On A j at t n +1 . Let α ( x, y, t n +1 ) and β ( x, y, t n +1 ) ∈ Q ( x, y ) interpolate a and b functions respectively at four vertices of A j , e.g. α ( x lj , y bj , t n +1 ) = a ( x lj , y bj , t n +1 ) , α ( x lj , y tj , t n +1 ) = a ( x lj , y tj , t n +1 ) , (3.3) α ( x rj , y bj , t n +1 ) = a ( x rj , y bj , t n +1 ) , α ( x rj , y tj , t n +1 ) = a ( x rj , y tj , t n +1 ) . Similarly, β is a bilinear function interpolating b at four vertices ( x lj , y bj ),( x lj , y tj ), ( x rj , y bj ), ( x rj , y tj ).0 X. CAI, J.-M. QIU, AND Y. YANG On ˜ A j ( t ) at t ∈ [ t n , t n +1 ) . Along characteristic lines of the adjoint problem(3.2) emanating from any point ( ξ, η ) ∈ A j at t n +1 , with˜ x ( t ; ( ξ, η, t n +1 )) , ˜ y ( t ; ( ξ, η, t n +1 ))satisfy the following equations,(3.4) ddt ˜ x ( t ; ( ξ, η, t n +1 )) = α ( ξ, η, t n +1 ) , ddt ˜ y ( t ; ( ξ, η, t n +1 )) = β ( ξ, η, t n +1 ) , from which one have(3.5) ˜ x ( τ ; ( ξ, η, t n +1 )) = ξ − α ( ξ, η, t n +1 )( t n +1 − τ ) ∈ Q ( ξ, η ) , (3.6) ˜ y ( τ ; ( ξ, η, t n +1 )) = η − β ( ξ, η, t n +1 )( t n +1 − τ ) ∈ Q ( ξ, η ) , with the Jacobian(3.7) J ( ξ, η, τ ) = ∂ (˜ x, ˜ y ) ∂ ( ξ, η ) ( τ ) = (cid:32) − ∂α∂ξ ( t n +1 − τ ) ∂α∂η ( t n +1 − τ ) − ∂β∂ξ ( t n +1 − τ ) 1 − ∂β∂η ( t n +1 − τ ) (cid:33) . Then we let, for t ∈ [ t n , t n +1 ], and (˜ x, ˜ y ) ∈ ˜ A j ( t ),(3.8) α (˜ x ( t ; ( ξ, η, t n +1 )) , ˜ y ( t ; ( ξ, η, t n +1 )) , t ) = α ( ξ, η, t n +1 ) , (3.9) β (˜ x ( t ; ( ξ, η, t n +1 )) , ˜ y ( t ; ( ξ, η, t n +1 )) , t ) = β ( ξ, η, t n +1 ) . It can be easily checked that, (3.4) are the characteristics equations for themodified adjoint problem (3.2) with α and β functions defined by eq. (3.8)and (3.9). For the modified adjoint problem, the solution ψ stays constantalong characteristics, therefore we have(3.10) ψ (˜ x ( τ ; ( ξ, η, t n +1 )) , ˜ y ( τ ; ( ξ, η, t n +1 )) , τ ) = Ψ( ξ, η ) ∈ P k ( A j ) , ∀ τ ∈ [ t n , t n +1 ] . Next we introduce a few notations and useful equalities [11, 26] regarding thecoordinate transformation defined by (3.5)-(3.6) .(3.11) dxdy = det( J ( ξ, η )) dξdη, (3.12) ∇ x,y ψ ( x, y ) = J ( ξ, η ) − T ∇ ξ,η Ψ( ξ, η ) , (3.13) n dS = det( J ( ξ, η )) J ( ξ, η ) − T ˘ n d ˘ S, where dS and d ˘ S are the infinitesimal boundaries of the dynamic element and theisoparametric element, respectively and their corresponding normal vectors are n and˘ n . The inverse of the Jacobian is given by(3.14) J ( ξ, η ) − = 1 | det( J ( ξ, η )) | (cid:18) ˜ y η − ˜ x η − ˜ y ξ ˜ x ξ (cid:19) . N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD ˜ A ( ) j ˜ A j ( t ) A j y xt ( x lj , y bj ) ( x rj , y bj )( x rj , y tj )( x lj , y tj ) η ξ Fig. 3 . Illustration for the mapping between dynamic element ˜ A j ( t ) (left) and the iso-parametric element (right). We assume the determinant of the Jacobian det( J ( ξ, η )) is positive; if the Jacobianis negative, it indicates the distortion of upstream cells. In such a situation, the timestepping size should be reduced by using the adaptive time stepping algorithm [4]. (2) Semi-discrete ELDG scheme formulation. Integrating (3.1) · ψ + (3.2) · u over Ω j , we have(3.15) (cid:90) Ω j [(3.1) · ψ + (3.2) · u ] dxdydt = 0 . Then,0 = (cid:90) t n +1 t n (cid:90) ˜ A j ( t ) ( u t ψ + uψ t ) dxdydt + (cid:90) t n +1 t n (cid:90) ˜ A j ( t ) (( au ) x ψ + αψ x u + ( bu ) y ψ + βψ y u ) dxdydt = (cid:90) t n +1 t n (cid:34)(cid:90) ˜ A j ( t ) ( uψ ) t dxdydt + (cid:90) ˜ A j ( t ) (( au ) x ψ + αψ x u + ( bu ) y ψ + βψ y u ) dxdydt (cid:35) = (cid:90) t n +1 t n (cid:34) ddt (cid:90) ˜ A j ( t ) uψdxdy − (cid:90) ∂ ˜ A j ( t ) uψ (cid:18) αβ (cid:19) · n dS (cid:35) dt + (cid:90) t n +1 t n (cid:34)(cid:90) ˜ A j ( t ) ∇ · (cid:18) aubu (cid:19) ψdxdy + (cid:90) ˜ A j ( t ) (cid:18) αβ (cid:19) · ∇ ψudxdy (cid:35) dt = (cid:90) t n +1 t n (cid:34) ddt (cid:90) ˜ A j ( t ) uψdxdy + (cid:90) ∂ ˜ A j ( t ) ψ F · n dS − (cid:90) ˜ A j ( t ) F · ∇ ψdxdy (cid:35) dt, (3.16)with(3.17) F ( u, x, y, t ) = (cid:18) ( a ( x, y, t ) − α ( x, y, t )) u ( b ( x, y, t ) − β ( x, y, t )) u (cid:19) , X. CAI, J.-M. QIU, AND Y. YANG in which the Leibniz-Reynolds transport theorem and the divergence Theorem areused for the above derivation. The time differential version of eq. (3.16) can bewritten as ddt (cid:90) ˜ A j ( t ) uψdxdy = − (cid:90) ∂ ˜ A j ( t ) ψ F · n dS + (cid:90) ˜ A j ( t ) F · ∇ ψdxdy. (3.18)As the 1D case, we map the coordinate of ( x, y ) ∈ ˜ A j ( t ) to a reference cell of( ξ, η ) ∈ A j as shown in Figure 3. Then we rewrite eq. (3.18) as ddt (cid:90) A j u (˜ x ( t, ( ξ, η, t n +1 )) , ˜ y ( t, ( ξ, η, t n +1 )) , t )Ψ( ξ, η ) det( J ( ξ, η, t )) dξdη = − (cid:90) ∂A j Ψ( ξ, η ) F · (cid:0) det( J ( ξ, η, t )) J ( ξ, η, t ) − T ˘ n (cid:1) d ˘ S + (cid:90) A j F · ( J ( ξ, η, t ) − T ∇ ξ,η Ψ) det( J ( ξ, η, t )) dξdη. (3.19)Notice that in equation (3.19), functions are all in the ( ξ, η ) coordinate, and can beevolved by the method-of-lines approach, e.g. using explicit SSP RK methods. Ψ( ξ, η )function stays as the same polynomial in the ( ξ, η ) coordinate for all t ∈ [ t n , t n +1 ] bythe design of our adjoint problem, see eq. (3.10).We let the approximate solution of u (˜ x ( t, ( ξ, η, t n +1 )) , ˜ y ( t, ( ξ, η, t n +1 )) , t ) be writ-ten in the ( ξ, η ) coordinate as follows,(3.20) u h ( ξ, η, t ) = n k (cid:88) p =1 ˘ u p ( t )Ψ p ( ξ, η ) , where bases Ψ p ( ξ, η ) , p = 1 , · · · , n k expands the space of P k ( A j ), for implementation.For the ELDG scheme, we look for u h in the above form satisfying ddt (cid:90) A j u h Ψ p det( J ( ξ, η, t )) dξdη + (cid:90) ∂A j Ψ p ˆ F · (cid:0) det( J ( ξ, η, t )) J ( ξ, η, t ) − T ˘ n (cid:1) d ˘ S − (cid:90) A j F · ( J ( ξ, η, t ) − T ∇ ξ,η Ψ) det( J ( ξ, η, t )) dξdη = 0 . (3.21)Here ˆ F in the second term is a monotone numerical flux, an example of which is theLax-Friedrich flux, and the line and volume integral in the second and third termscould be performed by proper high order quadrature rules as in a standard RK DGscheme. Then the coefficients u = (˘ u , ˘ u , · · · , ˘ u n k ) T in (3.20) satisfies a system ofODEs,(3.22) ddt ( M ( t ) u ( t )) = L ( u ( t )) , where the mass matrix M is of size n k by n k and its entries are M pq ( t ) = (cid:90) A j Ψ p ( ξ, η )Ψ q ( ξ, η ) det( J ( ξ, η, t )) dξdη, and L ( u ( t )) is the RHS vector from the evaluation of the other terms in (3.21). N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD (3) RK time discretization and fully discrete scheme. The semi-discretescheme (3.22) can be discretized by applying an explicit RK time discretization withthe initial condition(3.23) M ( t n ) u n = (cid:90) A (cid:63)j u nh ( x, y ) ψ ( x, y, t n ) dxdy, being evaluated by a 2D SLDG procedure [5]. Below we provide a flow chart of thefully discrete 2D algorithm described above. Step 1.
Construct α ( x, y, t ) and β ( x, y, t ) for ( x, y, t ) ∈ ˜ A j ( t ) × [ t n , t n +1 ] by firstconstructing α ( x, y, t n +1 ) , β ( x, y, t n +1 ) ∈ Q ( x, y ) , ( x, y ) ∈ A j , interpolating a ( x, y, t n +1 ) , b ( x, y, t n +1 ) respectively at four vertices of A j ;then these α and β functions are constructed by following (3.8)-(3.9) for t ∈ [ t n , t n +1 ). In particular, one first find ( ξ, η ) for (˜ x , ˜ y ) from (3.5)-(3.6);then the α (˜ x , ˜ y, t ) and β (˜ x , ˜ y, t ) are defined following (3.8)-(3.9). Note that,while (˜ x, ˜ y ) is a bilinear function of ( ξ, η ), the same statement does not holdfor the inverse mapping. Figure 3 illustrates 2D transformation between( ξ, η ) ∈ A j and ( x, y ) ∈ ˜ A j ( t ) for some t ∈ [ t n , t n +1 ]. Step 2.
Set up dynamic elements ˜ A ( l ) j , l = 0 , · · · , s, for each immediate stage of theRK method, and compute the corresponding Jacobian of the transformation J = ∂ ( x,y ) ∂ ( ξ,η ) , J ( ξ, η, τ ) − in (3.19); these quantities can be precomputed asfunctions of ( ξ, η, t ( l ) ). Step 3.
Perform the SLDG algorithm in [5] to get the initial condition of (3.23).Notice that since the mapping ( x ( ξ, η ) , y ( ξ, η )) in (3.5)-(3.6) is not affine, itis not as straighforward to find the inverse mapping of ( ξ ( x, y ) , η ( x, y )) as the1D problem. Some approximation, as is done in [5], has to be performed inorder to obtain ψ ( x, y, t n ). Step 4.
An SSP RK method is applied to (3.22). In particular, at the l th RK stage, M ( l ) u ( l ) is first being updated, then u ( l ) is computed by applying ( M ( l ) ) − ;finally u ( l ) as the degree of freedom in ( ξ, η ) coordinate are being used toevaluate the RHS of (3.22) for future RK stages. Remark α ( x, y, t n +1 )and β ( x, y, t n +1 ) functions are in Q ( A j ) in the modified adjoint problem ensures thequadrilateral shape of upstream cells. This avoids the need to use quadratic curvesto approximate upstream cells in achieving high order spatial accuracy in the originalSLDG algorithm [5]. An example of such is the swirling deformation example asshown in the numerical section. Remark a ( x, y, t ) , b ( x, y, t )) that are smoothenough and divergence free. The proposed ELDG formulation works for general non-divergence free velocity field as long as the Jacobian of the transformation is alwayspositive.
4. ELDG method with the exponential integrators for nonlinear Vlasovdynamics.
The proposed ELDG method for linear transport problems can be appliedto solve nonlinear models such as Vlasov models, via combining with the Runge-Kuttaexponential integrator method in [10, 4]. We will denote such a method as ELDG-RKEI. Below we first present the nonlinear Vlasov-Poisson, the guiding center Vlasov4
X. CAI, J.-M. QIU, AND Y. YANG models as well as the 2D incompressible Euler equations; and then present a secondorder and a third order ELDG-RKEI method.
The nonlinear Vlasov-Poisson system reads as follows,(4.1) f t + vf x + E ( x, t ) f v = 0 , (4.2) E ( x, t ) = − φ x , − φ xx ( x, t ) = ρ ( x, t ) , where the electron distribution function f ( x, v, t ) is the probability distribution func-tion in the phase space ( x, v ) ∈ Ω x × R describing the probability of finding a particlewith velocity v at position x and at time t . The electric field E = − φ x , where theself-consistent electrostatic potential φ is determined by the Poisson’s equation (4.2). ρ ( x, t ) = (cid:82) R f ( x, v, t ) dv − The guiding center Vlasov model describes a highly magnetized plasma in thetransverse plane of a tokamak [30, 16], and reads as follows:(4.3) ρ t + ∇ · ( E ⊥ ρ ) = 0 , (4.4) − ∆Φ = ρ, E ⊥ = ( − Φ y , Φ x ) , where the unknown variable ρ denotes the charge density of the plasma, and theelectric field E depends on ρ via the Poisson equation. The 2D incompressible Euler in the vorticity-stream function reads as fol-lows,(4.5) ω t + ∇ · ( u ω ) = 0 , (4.6) ∆Φ = ω, u = − ( − Φ y , Φ x ) , where u is the velocity field, ω is the vorticity of the fluid, and Φ is the stream-functiondetermined by Poisson’s equation.The above three models can be written in the form of (1.1). In [10, 9, 4], theexponential integrator method is applied to solve nonlinear time-dependent problems(1.1), by decomposing the nonlinear dynamics into the composition of a sequence oflinearized transport problems to achieve high order temporal accuracy. We denotethe ELDG procedure of updating the solution of linearized equation from t ∗ to t ∗ +∆ t with frozen velocity field P ( u ∗ ; x , t ∗ )(4.7) (cid:40) u t + ∇ · ( P ( u ∗ ; x , t ∗ ) u ) = 0 ,u ( t ∗ ) = u ∗ , as(4.8) ELDG ( P ( u ∗ ; x , t ∗ ) , ∆ t )( u ∗ ) . When a second order RKEI scheme is used with the ELDG update of linearizedsolution, one has u (1) = u n u (2) = ELDG (cid:18) P ( u (1) ) , ∆ t (cid:19) u (1) u n +1 = ELDG (cid:16) P ( u (2) ) , ∆ t (cid:17) u (1) . N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD u (1) = u n u (2) = ELDG (cid:18) P ( u (1) ) , ∆ t (cid:19) u (1) u (3) = ELDG (cid:18) P ( u (2) ) , ∆ t (cid:19) u (1) u n +1 = ELDG (cid:18) − P ( u (1) ) + 34 P ( u (3) ) , ∆ t (cid:19) u (2) . We name such scheme ‘ELDG-CF3C03’ [4], in which ‘CF3C03’ refers to the abovethird order RKEI scheme. We refer to [4] for more details regarding implementation.In the nonlinear Vlasov models LDG schemes [1, 14, 7, 28] are adopted to solve theelliptic field equations (4.2) and (4.4).
5. Numerical results.
In this section, we perform numerical experiments forlinear transport problems as well as the nonlinear Vlasov models. To showcase theproposed method, we perform the following studies: (1) the convergence of spatialdiscretization by using small enough time stepping size; (2) we vary
CF L to studythe temporal convergence and numerical stability with a well resolved spatial mesh;(3) we present snapshots of numerical solutions in a long time; (4) we numericallytrack the time history of invariants, such as mass and energy.The ELDG method presented below is the ELDG-ST1 method, unless otherwisenoted. When needed, we use the k + 1-th order RK for tracing characteristic lines.We set the time step for 1D and 2D problems as(5.1) ∆ t = CF L ∆ x and ∆ t = CF L a ∆ x + b ∆ y , respectively; here a and b are maximum transport speeds in x and y directions, re-spectively. For some test cases, we also present the SLDG [5, 4] and classical RKDGmethods for comparison purpose. Example (1D linear transport equation with constant coefficient.) We startwith the following 1D transport equation (5.2) u t + u x = 0 , x ∈ [0 , π ] , with the smooth initial data u ( x,
0) = sin( x ) and exact solution u ( x, t ) = sin( x − t ) .For the constant coefficient problem, the proposed ELDG method, if using the exactvelocity field, is the same as SLDG. Here we perturb the velocity at cell boundariesfor the modified adjoint problem to be α ( x j + ) = 1 + sin( x j + )∆ x .Table 2 reports the spatial accuracies of the ELDG, SLDG and RKDG methodsfor this example with the same time stepping size. The proposed ELDG method isfound to be as accurate as the SLDG and RKDG methods. We vary time steppingsize, with fixed well-resolved spatial meshes, and plot error vs. CF L in Figure 4 forELDG and SLDG P (left) and P (right) schemes at a long time T = 100 . For X. CAI, J.-M. QIU, AND Y. YANG the ELDG scheme, the time-stepping constraint can be found to be ∆ t ≤ k +1)∆ x ∆ x from the perturbation of velocity field and (2.22) ; hence CF L upper = 1(2 k + 1)∆ x , for P k ELDG schemes. They are shown as dashed lines in the figure. It is observedthat these bounds are expected in this numerical test. The SLDG schemes are observedto be unconditionally stable. The ELDG and SLDG schemes are observed to havesimilar error magnitudes, when the
CF L is less than the stability bounds (dash lines).
Table 2
1D linear transport equation with constant coefficient. u t + u x = 0 with initial condition u ( x,
0) = sin( x ) . T = π . We use CF L = 0 . and CF L = 0 . for all P and P schemes,respectively. ELDG here with the vertex perturbation. Mesh L error Order L error Order L error Order P RKDG P SLDG P ELDG40 1.15E-03 – 6.37E-04 – 6.08E-04 –80 2.85E-04 2.01 1.59E-04 2.00 1.55E-04 1.97160 7.09E-05 2.01 3.90E-05 2.03 3.84E-05 2.02320 1.77E-05 2.00 1.77E-05 2.00 9.77E-06 1.98 P RKDG P SLDG P ELDG40 9.28E-06 – 7.25E-06 – 7.69E-06 –80 1.16E-06 3.00 9.23E-07 2.97 9.45E-07 3.03160 1.44E-07 3.00 1.17E-07 2.98 1.18E-07 3.00320 1.80E-08 3.00 1.40E-08 3.06 1.41E-08 3.07
Fig. 4 . The L ∞ error versus CF L of SLDG methods and ELDG methods for 1D linear trans-port equation with constant coefficient: u t + u x = 0 with initial condition u ( x,
0) = sin( x ) . A longtime simulation is performed with T = 100 . The vertical long dashes from left to right are expectedupper bounds of CF L for stability for P k ELDG methods with meshes , and respectively. Example (1D transport equation with variable coefficients.) Consider (5.3) u t + (sin( x ) u ) x = 0 , x ∈ [0 , π ] with initial condition u ( x,
0) = 1 and the periodic boundary condition. The exactsolution is given by (5.4) u ( x, t ) = sin(2 tan − ( e − t tan( x )))sin( x ) . N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD As in the previous example, the spatial convergence of RKDG, SLDG, ELDG-ST1and ELDG-ST2 are shown in Table 3. The expected spatial convergence orders areobserved. In Figure 5, we plot the L ∞ error versus CF L of ELDG-ST1, ELDG-ST2and SLDG schemes with P (left) and P (right) polynomial spaces. The followingobservations are made: (1) all methods perform similarly around and before CF L = 1 ,which is well above the stability constraint of the RKDG method / (2 k + 1) ; (2) after CF L = 1 and before stability constraint of the method, the temporal convergence orderis observed to be and for P and P respectively, corresponding to the RK methodused in time integration and characteristics tracing; (3) the upper bounds of CF L forstability of P ELDG with mesh N = 80 , , are around . , , , which increasewith ratio around √ . This verifies the time step estimate ∆ t ∼ √ ∆ x in Remark 2.5. Table 3
1D transport equation with variable coefficients. u t + (sin( x ) u ) x = 0 with the initial condition u ( x,
0) = 1 . T = 1 . We use CF L = 0 . and CF L = 0 . for all P and P schemes, respectively. Mesh L error Order L error Order L error Order L error Order P RKDG P SLDG P ELDG-ST1 P ELDG-ST240 1.30E-03 – 1.35E-03 – 1.20E-03 – 1.35E-03 –80 3.25E-04 2.00 3.56E-04 1.92 3.24E-04 1.89 3.54E-04 1.93160 8.14E-05 2.00 8.95E-05 1.99 8.35E-05 1.96 8.89E-05 1.99320 2.04E-05 2.00 2.31E-05 1.95 2.21E-05 1.92 2.30E-05 1.95 P RKDG P SLDG P ELDG-ST1 P ELDG-ST240 8.11E-05 – 5.16E-05 – 6.45E-05 – 5.20E-05 –80 1.21E-05 2.74 6.35E-06 3.02 7.36E-06 3.13 6.36E-06 3.03160 1.79E-06 2.76 7.85E-07 3.02 8.65E-07 3.09 7.87E-07 3.02320 2.62E-07 2.78 9.61E-08 3.03 1.02E-07 3.08 9.63E-08 3.03
Fig. 5 . The L ∞ error versus CF L of SLDG methods and ELDG methods for 1D transportequation with variable coefficients. u t + (sin( x ) u ) x = 0 with the initial condition u ( x,
0) = 1 . T = 1 . ∆ t = CF L ∆ x . P SLDG-E means P SLDG scheme tracking characteristic lines exactly.
Example (Rigid body rotation.) Consider (5.5) u t − ( yu ) x + ( xu ) y = 0 , ( x, y ) ∈ [ − π, π ] . X. CAI, J.-M. QIU, AND Y. YANG
The initial condition is set to be the following smooth cosine bell (with C smoothness), (5.6) u ( x, y,
0) = (cid:40) r b cos (cid:16) r b r b π (cid:17) , if r b < r b , , otherwise , where r b = 0 . π , and r b = (cid:112) ( x − x b ) + ( y − y b ) denotes the distance between ( x, y ) and the center of the cosine bell ( x b , y b ) = (0 . π, . First of all, we present thespatial accuracies of ELDG, SLDG and RKDG for solving this problem up to T = 2 π in Table 4; the expected k + 1 th order of convergence is observed for these schemeswith P k polynomial space. Then, we study numerical stabilities of ELDG and SLDGmethods. In Figure 6, we present the plots of L ∞ error versus CF L of ELDG andSLDG schemes with different meshes. A few observations can be made: (1) When
CF L is around and below order , both schemes have similar performance in errormagnitude and order of convergence. Notice that this time stepping size is well abovethe stability constraint of / (2 k + 1) for RKDG. (2) When CF L is relatively largebut smaller than the stability constraint of ELDG, the temporal error starts to kickin 2nd and 3rd order temporal convergence order is shown. (3) Maximum
CF L s of P ELDG-ST1 using N = 40 , , are around , , . The increasing rate isaround . . Maximum CFLs of P ELDG-ST2 using N = 40 , , are around , . , . . The increasing rate is around . . The increasing ratio of upper boundsof CF L is around √ , which coincides with ∆ t ∼ √ ∆ x as in Remark 2.5. Similarobservations can be made for the P case. Table 4
Rigid body rotation. u t − ( yu ) x + ( xu ) y = 0 with the smooth cosine bell. T = 2 π. We use
CF L = 0 . and CF L = 0 . for all P and P schemes, respectively. Mesh L ∞ error Order L ∞ error Order L ∞ error Order L ∞ error Order P RKDG P SLDG P ELDG-ST1 P ELDG-ST220 P RKDG P SLDG-QC P ELDG-ST1 P ELDG-ST220 Example (Swirling deformation flow.) We consider solving (5.7) u t − (cid:16) cos (cid:16) x (cid:17) sin( y ) g ( t ) u (cid:17) x + (cid:16) sin( x ) cos (cid:16) y (cid:17) g ( t ) u (cid:17) y = 0 , ( x, y ) ∈ [ − π, π ] , with the same initial condition (5.6) , where g ( t ) = cos (cid:0) πtT (cid:1) π and T = 1 . . AsExample 5.3, we also study the spatial error and the numerical stability of the proposedELDG schemes in Table 5 and Figure 7, respectively. The similar observations asExample 5.3 can be made. Example (Vlasov-Poisson system: strong Landau damping.) Consider thestrong Landau damping for the Vlasov-Poisson system (4.1) with the initial condition
N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD Fig. 6 . The L ∞ error versus CF L of SLDG schemes and ELDG schemes for the rigid bodyrotation with the smooth cosine bells (5.6) . T = 2 π . Table 5
Swirling deformation flow with the smooth cosine bells (5.6) . T = 1 . . We use CF L = 0 . and CF L = 0 . for all P and P schemes, respectively. Mesh L ∞ error Order L ∞ error Order L ∞ error Order L ∞ error Order P RKDG P SLDG P ELDG-ST1 P ELDG-ST220 P RKDG P SLDG-QC P ELDG-ST1 P ELDG-ST220 Fig. 7 . The L ∞ error versus CF L of SLDG methods and ELDG methods for the swirlingdeformation flow with the smooth cosine bells (5.6) with T = 1 . . being a perturbed equilibrium (5.8) f ( x, v, t = 0) = 1 √ π (1 + α cos( kx )) exp (cid:18) − v (cid:19) , with α = 0 . and k = 0 . on a computational domain, [0 , π ] × [ − π, π ] . There X. CAI, J.-M. QIU, AND Y. YANG are several invariants of this problem which should remain constant in time. Theseinclude L p norms, kinetic energy and entropy: • L p norm, ≤ p ≤ ∞ : (5.9) (cid:107) f (cid:107) p = (cid:18)(cid:90) v (cid:90) x | f ( x, v, t ) | p dxdv (cid:19) p , • Energy: (5.10)
Energy = (cid:90) v (cid:90) x f ( x, v, t ) v dxdv + (cid:90) x E ( x, t ) dx, • Entropy: (5.11)
Entropy = (cid:90) v (cid:90) x f ( x, v, t ) log( f ( x, v, t )) dxdv. This is a classical problem that has been numerically investigated by several authors,e.g. see [33, 38, 22, 6].We first test the spatial accuracy of ELDG with the third order temporal schemefor this problem and report the results in Table 6. The time reversibility of the Vlasov-Poisson system [17] is used to test the order of convergence. In Table 6, we showthe L errors and the corresponding orders of convergence for P k ELDG and SLDG, k = 1 , with CF L = 0 . . We observe the expected orders of convergence of ELDGand SLDG.We then test the numerical stability of ELDG schemes with different meshes forthis problem integrated to T = 5 . Figure 8 reports L ∞ errors versus CF L of solutionsof ELDG schemes as well as the SLDG scheme. From this Figure, we find the expectedorders of convergence of the temporal schemes; we also find that the scheme can allowfor as large as
CF L = 50 ; we observe that the results of ELDG are very close to thoseof SLDG.We next study the performances of ELDG for conserving invariants of this prob-lem. The parameters of the tests are set as follows: we use a mesh of × cellsand CF L = 10 . For mass conservation, we observed that the mass deviation of ELDGschemes is around − × − due to the domain cut-off in the velocity space; we omitthis result. Figure 9 shows time evolutions of the relative deviation of L norms of thesolution as well as the discrete kinetic energy and entropy. We make the observationsfor this Figure: P ELDG performs better than P ELDG for conserving L norm,as SLDG schemes; for conserving energy, ELDG is worse than SLDG; for conservingentropy, ELDG does a better job than SLDG.Finally, we study ELDG schemes for this problem for a long-time simulation.We present the plots of solutions of ELDG schemes at T = 40 in the middle andright panels of Figure 8. We observe that P ELDG performs much better than P ELDG for capturing the filamentation structures. We find that the solutions of both P and P ELDG are negative around the places where the density is close to vacuum.Therefore, the positivity-preserving limiter should be added to the current scheme, forwhich we plan to explore in the future.
Example (The guiding center Vlasov model: spatial accuracy and conver-gence test.) Consider the guiding center Vlasov model on the domain [0 , π ] × [0 , π ] N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD Table 6
Strong Landau damping. T = 0 . . Use the time reversibility of the VP system. Order ofaccuracy in space for the SLDG method and the ELDG method. The third order temporal schemeCF3C03 is used for all schemes. We set CF L = 0 . so that the spatial error is the dominant error. Mesh L error Order L error Order L error Order L error Order P SLDG P ELDG P SLDG-QC P ELDG32 Fig. 8 . Left panel: plots of L ∞ errors versus the CF L number for solving Strong Landaudamping at T = 5 . Temporal order of convergence in L ∞ norm of ELDG schemes as well asthe SLDG scheme coupled with exponential integrators by comparing numerical solutions with areference solution from the corresponding scheme with CF L = 0 . .Middle and right panels: surface plots of the numerical solutions for the strong Landau damping at T = 40 . We use a mesh of × cells and CF L = 10 . Middle: P ELDG+CF2. Right: P ELDG+CF3C03.
Fig. 9 . Strong Landau damping. Time evolutions of the relative deviation of L (left) normsof the solution as well as the discrete kinetic energy (middle) and entropy (right). We use a meshof × cells and CF L = 10 for all simulations. with the initial condition, ρ ( x, y,
0) = − x ) sin( y ) and the periodic boundary con-dition. The exact solution stays stationary. We test the spatial convergence of the pro-posed ELDG schemes as well as SLDG schemes with the third order temporal scheme,CF3C03, for solving the guiding center Vlasov model up to time T = 1 and reportthe results in Table 7. We make the following observations: (1) we find the expectedorders of convergence for P k ELDG+ P k +1 LDG, k = 1 , , in L and L ∞ norms; (2)the results of ELDG schemes are almost the same as those of SLDG schemes. Example (The guiding center Vlasov model: Kelvin-Helmholtz instabilityproblem.) We consider the two-dimensional guiding center model problem (4.3) with X. CAI, J.-M. QIU, AND Y. YANG
Table 7
The guiding center Vlasov model on the domain [0 , π ] × [0 , π ] with the initial condition ρ ( x, y,
0) = − x ) sin( y ) . Periodic boundary conditions in two directions. Spatial orders of con-vergence of P k SLDG(-QC)+ P r LDG+CF3C03 and P k ELDG+ P r LDG+CF3C03, k = 1 , , and r = k + 1 . T = 1 . CF L = 1 . Mesh L error Order L ∞ error Order L error Order L ∞ error Order P SLDG P ELDG20 P SLDG-QC P ELDG20 the initial condition (5.12) ρ ( x, y ) = sin( y ) + 0 .
015 cos( kx ) , and periodic boundary condition on the domain [0 , π ] × [0 , π ] . We let k = 0 . , whichwill create a Kelvin-Helmholtz instability [30].First, we test the temporal convergence of the proposed ELDG schemes with dif-ferent temporal schemes by computing this problem up to T = 5 . In particular, we testthe proposed second scheme, P ELDG+ P LDG+CF2, and the third order scheme, P ELDG+ P LDG+CF3C03. In order to minimize the errors for the spatial scheme,a fixed mesh of × cells is used. The reference solution is computed by the samescheme with the same mesh but using a small CF L = 0 . . We show the plots of L er-rors versus the CF L number of the proposed ELDG schemes for the Kelvin-Helmholtzinstability problem at T = 5 in Figure 10. We make a few observations: (1) we ob-serve expected orders of convergence for all temporal schemes; and CF L of ELDG canbe taken to be as large as 50; (2) by comparing the error magnitude, P ELDG+ P LDG+CF3C03 performs slightly better than P SLDG-QC+ P LDG+CF3C03.
Fig. 10 . Plots of L errors versus the CF L number of the proposed ELDG schemes as wellas the SLDG scheme for the Kelvin-Helmholtz instability problem at T = 5 . Temporal order ofconvergence of presented schemes by comparing numerical solutions with a reference solution fromthe corresponding scheme with CF L = 0 . . The mesh of × cells is used. N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD We then study the quality of the proposed ELDG schemes by tracking relativedeviations of some invariants of this problem, the energy (cid:107) E (cid:107) L = (cid:82) Ω E · E dxdy andthe enstrophy (cid:107) ρ (cid:107) L = (cid:82) Ω ρ dxdy . We study ELDG schemes using a mesh of × cells with CF L = 5 for solving this problem for a long-time simulation and report theresults in Figure 11. We find that P ELDG can perform much better than P ELDGfor conserving both energy and enstrophy. We find that by comparing SLDG andELDG with the same polynomial space for conserving both energy and enstrophy, thecomparable results can be observed. Finally, we show surface plots of the numericalsolutions for the Kelvin-Helmholtz instability at T = 40 in Figure 12. We still observethat the resolution of solutions of ELDG is comparable to that of SLDG. Fig. 11 . Time evolutions of the relative deviation of energy (left) and enstrophy (right) for theproposed ELDG schemes for the Kelvin-Helmholtz instability problem. The mesh of × cellsand CF L = 5 are used.
Fig. 12 . Surface plots of the numerical solutions for the Kelvin-Helmholtz instability at T = 40 .We use a mesh of × cells and CF L = 5 . Left: P SLDG-QC+ P LDG+CF3C03. Right: P ELDG+ P LDG+CF3C03.
Example (The incompressible Euler equations: the shear flow problem) Forthe double shear layer problem [2, 36], we solve the 2D incompressible Euler equations (4.5) in the domain [0 , π ] × [0 , π ] , with the periodic boundary conditions and theinitial condition given by (5.13) ω ( x, y,
0) = δ cos( x ) − ρ sech (cid:16) y − π/ ρ (cid:17) , if y ≤ π,δ cos( x ) + ρ sech (cid:16) π/ − yρ (cid:17) , if y > π, where δ = 0 . and ρ = π/ . X. CAI, J.-M. QIU, AND Y. YANG
As time evolves, the solution quickly rolls up with smaller and smaller spatialscales so on any fixed grid, the full resolution will be lost eventually. This problemis a classic benchmark for demonstrating the effectiveness of a new scheme so it hasbeen tested for many schemes such as the high order nonsplitting SL WENO scheme[34], the DG method in [24, 36, 39] and the spectral element method in [18, 35]. Wefirst show surface plots of numerical solutions for this problem at T = 8 in Figure 13,where the solution is rolled up in a very small scale. We find that ELDG schemescould allow for CF L = 5 for these simulations and the solutions with larger
CF L = 5 seem to be less dissipative than those with
CF L = 1 . We then study the quality of theELDG schemes by tracking relative deviations of the energy (cid:107) u (cid:107) L = (cid:82) Ω u · u dxdy andthe enstrophy (cid:107) ω (cid:107) L = (cid:82) Ω ω dxdy of this problem and report the results in Figure 14.We observed that higher order P ELDG performs much better than the lower order P ELDG for conserving energy and enstrophy.
Fig. 13 . Contour plots of the numerical solutions for the shear flow test at T = 8 . P ELDG+ P LDG+CF3C03 using
CF L = 1 (left),
CF L = 5 (right). The mesh of × . Fig. 14 . Time evolution of the relative deviation of energy and enstrophy for the proposedELDG schemes for the shear flow test. Left: ELDG+ P LDG+CF2. Right: P ELDG + P LDG+CF3C03. We use a mesh of × and CF L = 5 .
6. Conclusion.
In this paper, we have developed a new Eulerian-Lagrangiandiscontinuous Galerkin (DG) method for transport problems. The new framework en-compasses the semi-Lagrangian DG and Eulerian Runge-Kutta DG in special cases;thus inherits advantages from both approaches in stability under large time step-ping sizes, and in mass conservation, compactness and high order accuracy. Theseadvantages are numerically verified by extensive numerical tests for linear transportequation and nonlinear dynamics. Future works include further theoretic development
N EULERIAN-LAGRANGIAN DISCONTINUOUS GALERKIN METHOD
REFERENCES[1]
D. Arnold, F. Brezzi, B. Cockburn, and L. Marini , Unified analysis of discontinuousGalerkin methods for elliptic problems , SIAM Journal on Numerical Analysis, 39 (2002),pp. 1749–1779.[2]
J. Bell, P. Colella, and H. Glaz , A second-order projection method for the incompressibleNavier-Stokes equations , Journal of Computational Physics, 85 (1989), pp. 257–283.[3]
P. A. Bosler, A. M. Bradley, and M. A. Taylor , Conservative Multimoment Transportalong Characteristics for Discontinuous Galerkin Methods , SIAM Journal on ScientificComputing, 41 (2019), pp. B870–B902.[4]
X. Cai, S. Boscarino, and J.-M. Qiu , High Order Semi-Lagrangian Discontinuous GalerkinMethod Coupled with Runge-Kutta Exponential Integrators for Nonlinear Vlasov Dynam-ics , arXiv preprint arXiv:1911.12229, (2019).[5]
X. Cai, W. Guo, and J.-M. Qiu , A high order conservative semi-Lagrangian discontinuousGalerkin method for two-dimensional transport simulations , Journal of Scientific Comput-ing, 73 (2017), pp. 514–542.[6]
X. Cai, W. Guo, and J.-M. Qiu , A high order semi-Lagrangian discontinuous Galerkinmethod for Vlasov-Poisson simulations without operator splitting , Journal of Computa-tional Physics, 354 (2018), pp. 529–551.[7]
P. Castillo, B. Cockburn, I. Perugia, and D. Sch¨otzau , An a priori error analysis ofthe local discontinuous Galerkin method for elliptic problems , SIAM Journal on NumericalAnalysis, 38 (2000), pp. 1676–1706.[8]
M. Celia, T. Russell, I. Herrera, and R. Ewing , An Eulerian-Lagrangian localized adjointmethod for the advection-diffusion equation , Advances in Water Resources, 13 (1990),pp. 187–206.[9]
E. Celledoni and B. K. Kometa , Semi-Lagrangian Runge-Kutta exponential integrators forconvection dominated problems , Journal of Scientific Computing, 41 (2009), pp. 139–164.[10]
E. Celledoni, A. Marthinsen, and B. Owren , Commutator-free Lie group methods , FutureGeneration Computer Systems, 19 (2003), pp. 341–352.[11]
P. G. Ciarlet , Mathematical Elasticity: Volume I: three-dimensional elasticity , North-Holland, 1988.[12]
B. Cockburn and C.-W. Shu , TVB Runge-Kutta local projection discontinuous Galerkinfinite element method for conservation laws II: general framework , Mathematics of Com-pututation, (1989), pp. 411–435.[13]
B. Cockburn and C.-W. Shu , The Runge-Kutta local projection-discontinuous-Galerkin fi-nite element method for scalar conservation laws , ESAIM: Mathematical Modelling andNumerical Analysis, 25 (1991), pp. 337–361.[14]
B. Cockburn and C.-W. Shu , The local discontinuous Galerkin method for time-dependentconvection-diffusion systems , SIAM Journal on Numerical Analysis, 35 (1998), pp. 2440–2463.[15]
B. Cockburn and C.-W. Shu , Runge–Kutta discontinuous Galerkin methods for convection-dominated problems , Journal of Scientific Computing, 16 (2001), pp. 173–261.[16]
N. Crouseilles, M. Mehrenberger, and E. Sonnendr¨ucker , Conservative semi-Lagrangianschemes for Vlasov equations , Journal of Computational Physics, 229 (2010), pp. 1927–1953.[17]
P. Degond, L. Pareschi, and G. Russo , Modeling and Computational Methods for KineticEquations , Springer, 2004.[18]
P. Fischer and J. Mullen , Filter-based stabilization of spectral element methods , ComptesRendus de l’Acad´emie des Sciences-Series I-Mathematics, 332 (2001), pp. 265–270.[19]
W. Guo, R. Nair, and J.-M. Qiu , A conservative semi-Lagrangian discontinuous Galerkinscheme on the cubed-sphere , Monthly Weather Review, 142 (2013), pp. 457–475.[20]
C.-S. Huang and T. Arbogast , An Eulerian–Lagrangian Weighted Essentially Nonoscilla-tory scheme for Nonlinear Conservation Laws , Numerical Methods for Partial DifferentialEquations, 33 (2017), pp. 651–680.[21]
C.-S. Huang and T. Arbogast , An Implicit Eulerian–Lagrangian WENO3 Scheme for Non-linear Conservation Laws , Journal of Scientific Computing, 77 (2018), pp. 1084–1114.[22]
C.-S. Huang, T. Arbogast, and C.-H. Hung , A semi-Lagrangian finite difference WENOscheme for scalar nonlinear conservation laws , Journal of Computational Physics, 322 X. CAI, J.-M. QIU, AND Y. YANG(2016), pp. 559–585.[23]
C. Klingenberg, G. Schn¨ucke, and Y. Xia , Arbitrary Lagrangian-Eulerian discontinuousGalerkin method for conservation laws: analysis and application in one dimension , Math-ematics of Computation, 86 (2017), pp. 1203–1232.[24]
J.-G. Liu and C.-W. Shu , A high-order discontinuous Galerkin method for 2D incompressibleflows , Journal of Computational Physics, 160 (2000), pp. 577–596.[25]
D. Luo, W. Huang, and J. Qiu , A quasi-Lagrangian moving mesh discontinuous Galerkinmethod for hyperbolic conservation laws , Journal of Computational Physics, 396 (2019),pp. 544–578.[26]
P.-O. Persson, J. Bonet, and J. Peraire , Discontinuous Galerkin solution of the Navier–Stokes equations on deformable domains , Computer Methods in Applied Mechanics andEngineering, 198 (2009), pp. 1585–1595.[27]
J.-M. Qiu and C.-W. Shu , Positivity preserving semi-Lagrangian discontinuous Galerkin for-mulation: Theoretical analysis and application to the Vlasov–Poisson system , Journal ofComputational Physics, 230 (2011), pp. 8386–8409.[28]
J. Rossmanith and D. Seal , A positivity-preserving high-order semi-Lagrangian discontinuousGalerkin scheme for the Vlasov-Poisson equations , Journal of Computational Physics, 230(2011), pp. 6203–6232.[29]
T. F. Russell and M. A. Celia , An overview of research on Eulerian–Lagrangian localizedadjoint methods (ELLAM) , Advances in Water resources, 25 (2002), pp. 1215–1231.[30]
M. M. Shoucri , A two-level implicit scheme for the numerical solution of the linearized vor-ticity equation , International Journal for Numerical Methods in Engineering, 17 (1981),pp. 1525–1538.[31]
C.-W. Shu and S. Osher , Efficient implementation of essentially non-oscillatory shock-capturing schemes , Journal of Computational Physics, 77 (1988), pp. 439–471.[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. Xiong, J.-M. Qiu, Z. Xu, and A. Christlieb , High order maximum principle preserv-ing semi-Lagrangian finite difference WENO schemes for the Vlasov equation , Journal ofComputational Physics, 273 (2014), pp. 618–639.[34]
T. Xiong, G. Russo, and J.-M. Qiu , High order multi-dimensional characteristics tracingfor the incompressible euler equation and the guiding-center vlasov equation , Journal ofScientific Computing, 77 (2018), pp. 263–282.[35]
C. Xu , Stabilization methods for spectral element computations of incompressible flows , Journalof Scientific Computing, 27 (2006), pp. 495–505.[36]
X. Zhang and C.-W. Shu , On maximum-principle-satisfying high order schemes for scalarconservation laws , Journal of Computational Physics, 229 (2010), pp. 3091–3120.[37]
L. Zhou, Y. Xia, and C.-W. Shu , Stability analysis and error estimates of arbitraryLagrangian-Eulerian discontinuous Galerkin method coupled with Runge-Kutta time-marching for linear conservation laws , ESAIM: Mathematical Modelling and NumericalAnalysis, 53 (2019), pp. 105–144.[38]
H. Zhu, J. Qiu, and J.-M. Qiu , An h-adaptive RKDG method for the Vlasov–Poisson system ,Journal of Scientific Computing, 69 (2016), pp. 1346–1365.[39]