A Generalized Eulerian-Lagrangian Discontinuous Galerkin Method for Transport Problems
AA Generalized Eulerian-Lagrangian Discontinuous Galerkin Methodfor Transport Problems ∗ Xue Hong † Jing-Mei Qiu ‡ February 24, 2021
Abstract
We propose a generalized Eulerian-Lagrangian (GEL) discontinuous Galerkin (DG) method.The method is a generalization of the Eulerian-Lagrangian (EL) DG method for transportproblems proposed in [arXiv preprint arXiv: 2002.02930 (2020)], which tracks solution alongapproximations to characteristics in the DG framework, allowing extra large time stepping sizewith stability. The newly proposed GEL DG method in this paper is motivated for solving linearhyperbolic systems with variable coefficients, where the velocity field for adjoint problems of thetest functions is frozen to constant. In this paper, in a simplified scalar setting, we propose theGEL DG methodology by freezing the velocity field of adjoint problems, and by formulating thesemi-discrete scheme over the space-time region partitioned by linear lines approximating char-acteristics. The fully-discrete schemes are obtained by method-of-lines Runge-Kutta methods.We further design flux limiters for the schemes to satisfy the discrete geometric conservationlaw (DGCL) and maximum principle preserving (MPP) properties. Numerical results on 1Dand 2D linear transport problems are presented to demonstrate great properties of the GEL DGmethod. These include the high order spatial and temporal accuracy, stability with extra largetime stepping size, and satisfaction of DGCL and MPP properties.
Key words:
Eulerian-Lagrangian; discontinuous Galerkin; characteristics method; mass conser-vative; discrete geometric conservation law; maximum principle preserving.
In this paper, we propose a generalized Eulerian-Lagrangian (GEL) Runge-Kutta (RK) discontin-uous Galerkin (DG) method for a model transport equation in the form of u t + ∇ · ( P ( u ; x , t ) u ) = 0 , ( x , t ) ∈ R d × [0 , T ] , (1.1) ∗ Research of the first author is supported by the China Scholarship Council for 2 years’ study at the University ofDelaware. Research of the second author is supported by NSF grant NSF-DMS-1818924, Air Force Office of ScientificResearch FA9550-18-1-0257 and University of Delaware. † School of Mathematical Sciences, University of Science and Technology of China, Hefei, Anhui, 230026, P.R.China. ([email protected]). ‡ Department of Mathematical Sciences, University of Delaware, Newark, DE, 19716, USA. ([email protected]). a r X i v : . [ m a t h . NA ] F e b here d is the spatial dimension, u : R d × [0 , T ] → R , and P ( u ; x , t ) = ( P ( u ; x , t ) , · · · , P d ( u ; x , t )) T is a linear or nonlinear velocity field. Such a model could come from a wide range of applicationfields including fluid dynamics, climate modeling, and kinetic description of plasma.The GEL DG method is a generalization from the EL DG method proposed in [5]. With RKtime discretization, their fully-discrete schemes are termed “GEL RK DG” and “EL RK DG”methods. The EL DG method is built upon a fixed set of computational mesh, yet in each timestep, the solution is evolved over a local dynamic space-time region Ω j (see Figure 2.1), the partitionof which is determined by linear approximations to characteristics. The EL DG method introducesmodified adjoint problems for the test function φ t + ˜P · ∇ φ = 0 , (1.2)with its velocity field ˜P being a linear approximation to the velocity field P in (1.1). Then RKmethods are used for the time discretization via the method-of-lines approach. The proposedGEL DG method shares the same space-time partition strategy as the EL DG method. A majordifference is the velocity field for the local modified adjoint problem of test functions. The GEL DGuses a constant function with ˜P ≡ ¯ P j in (1.2), whereas ˜P ∈ P ( x, t ) in the EL DG, to approximatethe velocity field P of (1.1). Such design is motivated from solving the wave equation via trackinginformation along different characteristics families in a linear system. Take a 1-D 2-by-2 hyperbolicsystem for example, U t + ( A ( x ) U ) x = 0 , which could arise from the wave propagation in a heterogeneous media. The modified adjointproblem for the system by the GEL DG method isΦ t + A j Φ x = 0 , where A j is a frozen local constant matrix that approximates A ( x ) on Ω j . In this paper, we focus onthe GEL RK DG algorithm for scalar transport problems. The proposed GEL RK DG maintainsmass conservation, high order spatial and temporal accuracy, and allows for extra large time stepswith stability. We also establish that the semi-discrete GEL DG and EL DG formulation aremathematically equivalent, whereas the time discretization introduces differences for fully discreteschemes. 2e further study the property of discrete geometric conservation law (DGCL) and maximumprinciple preserving (MPP) property of GEL RK DG method, and find that the method fail tosatisfy the DGCL and MPP property in general. We then propose MPP limiters to preserve theDGCL and MPP properties. The MPP limiters involve the polynomial rescaling limiter [26, 27]and the parametrized MPP flux limiter [22, 24, 23]. The polynomial rescaling limiter preservesthe MPP property for the piecewise DG polynomials, while the parametrized MPP flux limiterpreserves the MPP property of cell averages in the final RK stage only, to avoid order reduction ofRK methods if limiters are applied to intermediate RK solutions.Finally, among different classes of EL methods in the literature, we would like to mention a fewclosely related ones. Eulerian-Lagrangian finite volume methods were introduced in [13] to handlenonlinearity for characteristic methods in the finite volume framework. The Eulerian LagrangianLocalized Adjoint Methods (ELLAM) [6] introduces an adjoint problem for the test function in thecontinuous finite element framework and has been applied to different problems [20, 18]. Comparedwith ELLAM, the EL DG, SL DG and EL RK DG [21, 1, 5] are being developed in the discontinuousGalerkin finite element framework. Another line of development, that is closely related to thiswork, is the Arbitrary Lagrangian Eulerian (ALE) DG method [15, 12]. Both EL DG and ALEDG evolve the DG solution on a dynamic moving mesh. The dynamic mesh movement of the ELDG approximates characteristics for the potential of using larger time stepping sizes with stability,whereas the mesh movement of ALE DG could come from tracking moving computational domainand/or better shock resolution. The formulation of EL DG comes from the introduction of a localmodified adjoint problem, whereas the ALE DG method is formulated through the coordinatetransform of test function on a reference domain.This paper is organized as follows. In Section 2, we develop the GEL DG for one-dimensional(1D) linear transport problems. We discuss numerical treatment of inflow boundary conditions andextensions to 2D problems by dimensional splitting. In Section 3, we establish the equivalence ofthe GEL DG and EL DG method in semi-discrete form. In Section 4, we study the DGCL andMPP properties of the GEL RK DG method and propose a MPP limiter to preserve DGCL andMPP for fully discrete schemes. In Section 5, the performance of the proposed method is shownthrough extensive numerical tests. Finally, concluding remarks are made in Section 6.3 GEL DG formulation for linear transport problems
We propose the GEL DG method, which differs from the EL DG method [5] in the design of themodified adjoint problem for test functions. In the EL DG method, the adjoint problem is uniquelydetermined from the partition of the space-time region Ω j ; while in the GEL DG method, theadjoint problem is independent from the partition of the space-time region. Such design of adjointproblems offers flexibility in handling hyperbolic systems when characteristic decomposition variesin space. We start from a 1D linear transport equation in the following form u t + ( a ( x, t ) u ) x = 0 , x ∈ [ x a , x b ] . (2.1)For simplicity, we assume periodic boundary conditions, and the velocity field a ( x, t ) is a continuousfunction of space and time. We perform a partition of the computational domain x a = x < x < · · · < x N + = x b . Let I j = [ x j − , x j + ] denote an element of length ∆ x j = x j + − x j − and define∆ x = max j ∆ x j . We define the finite dimensional approximation space, V kh = { v h : v h | I j ∈ P k ( I j ) } ,where P k ( I j ) denotes the set of polynomials of degree at most k . We let t n be the n -th time leveland ∆ t = t n +1 − t n to be the time-stepping size.The scheme formulation is summarized in four steps. We will first partition the space-timeregion Ω j ’s, then introduce a new modified adjoint problem for the test function ψ ( x, t ). We thenformulate the semi-discrete GEL DG scheme. Finally, we apply the method-of-lines RK methodfor the time marching. (1) Partition of space-time region Ω j : We define a space-time region Ω j = ˜ I j ( t ) × [ t n , t n +1 ]with ˜ I j ( t ) = [˜ x j − ( t ) , ˜ x j + ( t )] , t ∈ [ t n , t n +1 ]being the dynamic interval, see Figure 2.1. Here ˜ x j ± ( t ) = x j ± + ( t − t n +1 ) ν j ± are straight linesemanating from cell boundaries x j ± with slopes ν j ± = a ( x j ± , t n +1 ) as in the EL DG scheme.We let I (cid:63)j . = ˜ I j ( t n ) = [ x ∗ j − , x ∗ j + ] be the upstream cell of I j at t n . Note that ν j ± are chosento best take advantage of the characteristics information. In a classical Eulerian RK DG scheme, ν j + = 0, ∀ j . 4 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 I ∗∗ j (cid:101) I j ( t ) ν j − ν j + α j α j ξ = x j − / ξ = x j +1 / Figure 2.1: Illustration for dynamic element (cid:101) I j ( t ) of new GEL DG. (2) Adjoint problems. We consider a local adjoint problem for the test function: ® ψ t + α j ψ x = 0 , ( x, t ) ∈ Ω j ,ψ ( t = t n +1 ) = Ψ( x ) , ∀ Ψ( x ) ∈ P k ( I ∗∗ j ) . (2.2)Here we let α j = a ( x j , t n +1 ). To obtain the test function ψ ( x, t ) on Ω j , Ψ( x ) needs to be definedin a large enough neighborhood containing I j , named I ∗∗ j = [ x ∗∗ j − , x ∗∗ j + ], with x ∗∗ j − = min( x j − , x ∗ j − + α j ∆ t ) , x ∗∗ j + = max( x j + , x ∗ j + + α j ∆ t ) . (2.3)Please see the green curves in Figure 2.1 for the slope of α j , and the interval I ∗∗ j . Here we take anatural extension of Ψ from I j to I ∗∗ j . The idea of defining Ψ( x ) on I ∗∗ j is to ensure that ψ ( x, t )can be properly found on Ω j , see the area shadowed by green dash lines in Figure 2.1. This isdifferent from the adjoint problem in the EL DG method, where the test function ψ stays the samepolynomial, if ˜ I j ( t ) is mapped to a reference interval I j . (3) Formulation of the semi-discrete GEL DG scheme. In order to formulate the scheme,we integrate (2.1) · ψ + (2.2) · u over Ω j , (cid:90) Ω j [(2.1) · ψ + (2.2) · u ] dxdt = 0 . (2.4)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 ψ + α j ψ 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 ) (( a ( x, t ) uψ ) x − a ( x, t ) uψ x + α j ψ x u ) dxdt = (cid:90) t n +1 t n ñ ddt (cid:90) ˜ I j ( t ) uψdx − ν j + uψ | ˜ x j + 12 ( t ) + ν j − uψ | ˜ x j − ( t ) + auψ (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t )˜ x j − ( t ) (cid:90) ˜ I j ( t ) ( α j − a ) uψ x dx ô dt = (cid:90) t n +1 t n ñ ddt (cid:90) ˜ I j ( t ) uψdx + ( a j + − ν j + ) uψ | ˜ x j + 12 ( t ) − ( a j − − ν j − ) uψ | ˜ x j − ( t ) + (cid:90) ˜ I j ( t ) ( α j − a ) uψ x dx ô dt. (2.5)Letting F ( u ) . = ( a − ν ) u , the time differential form of (2.5) gives ddt (cid:90) ˜ I j ( t ) ( uψ ) dx = − ( F ψ ) (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + ( F ψ ) (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t ) + (cid:90) ˜ I j ( t ) ( a − α j ) uψ x dx. (2.6)Notice that the dynamic interval of ˜ I j ( t ) can always be linearly mapped to a reference cell ξ in I j by the mapping ˜ x ( t ; ( ξ, t n +1 )) . then eq. (2.6) in the ξ -coordinate becomes ddt (cid:90) I j ( uψ ( ξ )) ∂ ˜ x ( t ; ( ξ, t n +1 )) ∂ξ dξ = − ( F ψ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ = x j + 12 + ( F ψ ) (cid:12)(cid:12)(cid:12)(cid:12) ξ = x j − + (cid:90) I j ( a − α j ) uψ ξ dξ. (2.7)The DG discretization of (2.7) is to find u h ( ξ, t ) ∈ P k ( I j ), so that the following equality holds, ddt (cid:90) I j u h ( ξ, t ) ψ ( ξ, t ) ∂ ˜ x∂ξ dξ = − ˆ F j + ψ (˜ x − j + ( t ) , t ) + ˆ F j − ψ (˜ x + j − ( t ) , t ) + (cid:90) I j ( a − α j ) u h ψ ξ dξ, (2.8)for ψ ( x, t ) satisfying the adjoint problem (2.2) with ∀ Ψ( x ) = ψ ( x, t n +1 ) ∈ P k ( I ∗∗ j ). Here ˆ F at acell boundary can be taken as a monotone flux, e.g. the Lax-Friedrichs fluxˆ F ( u − , u + ) = 12 ( F ( u − ) + F ( u + )) + α u − − u + ) , α = max u | F (cid:48) ( u ) | ; (2.9)and we use Gauss quadrature rules with k + 1 quadrature points to approximate the integral termon the R.S.H. of the equation (2.8).Next we discuss the choice of basis functions for representing solutions and test functions, withwhich one can assemble time-dependent mass matrices for implementation. In a classical EulerianDG setting, the test functions are the same as basis functions of V kh . However, in the GEL DGsetting, the test function ψ j,m ( x, t ) is in a time-dependent domain ˜ I j satisfying the adjoint problem(2.2) with Ψ( x ) = Ψ j,m ( x ) , j = 1 , ..., N, m = 0 , ..., k. (2.10) { Ψ j,m ( x ) } ≤ j ≤ N, ≤ m ≤ k are the basis of P k ( I j ) with a natural extension to I ∗∗ j . In fact, we havefrom eq. (2.2) ψ j,m ( x, t ) = Ψ j,m ( x − α j ( t − t n +1 )) . (2.11)6e let ˜ U j ( t ) be a vector of size ( k + 1) × ® (cid:90) ˜ I j ( t ) u h ( x, t ) ψ j,m ( x, t ) dx ´ ≤ m ≤ k . (2.12)On the other hand, we set our basis function as ¶ ˜ ψ j,m ( x, t ) © ≤ j ≤ N, ≤ m ≤ k on a reference cell witha mapping ˜ x ( t ;( ξ,t n +1 )) − ˜ x j − ( t )˜ x j + 12 ( t ) − ˜ x j − ( t ) = ξ − x j − ∆ x j with˜ ψ j,m (˜ x ( t ; ( ξ, t n +1 )) , t ) = Ψ j,m ( ξ ) , (2.13)where Ψ j,m | I j is a set of basis in V kh on I j . Then we let u h ( x, t ) = k (cid:88) l =0 ˆ u ( l ) j ( t ) ˜ ψ j,l ( x, t ) , on ˜ I j ( t ) , (2.14)where ˆ u ( l ) are coefficients for the basis. Let U j ( t ) = (ˆ u (0) j ( t ) , · · · , ˆ u ( k ) j ( t )) T be the coefficient vectorof size ( k + 1) ×
1. Notice that U j ( t ) here is different from ˜ U j ( t ) defined in (2.12), satisfying theGEL DG scheme (2.8) (cid:90) ˜ I j ( t ) u h ( x, t ) ψ j,m ( x, t ) dx = (cid:90) ˜ I j ( t ) k (cid:88) l =0 ˆ u ( l ) j ( t ) ˜ ψ j,l ( x, t ) ψ j,m ( x, t ) dx = k (cid:88) l =0 ˆ u ( l ) j ( t ) (cid:90) ˜ I j ( t ) ˜ ψ j,l ( x, t ) ψ j,m ( x, t ) dx = k (cid:88) l =0 ˆ u ( l ) j ( t ) (cid:90) ˜ I j ( t ) ˜ ψ j,l ( x, t )Ψ j,m ( x − α j ( t − t n +1 )) dx. (2.15)Now we assemble the time dependent mass-matrix A j ( t ) of size ( k + 1) × ( k + 1) with its elementsconsisting of ® (cid:90) ˜ I j ( t ) ˜ ψ j,l ( x, t )Ψ j,m ( x − α j ( t − t n +1 )) dx ´ ≤ l ≤ k, ≤ m ≤ k . Further, we have by eq. (2.12), (2.14), (2.15)˜ U j ( t ) = A j ( t ) · U j ( t ) , ∀ j = 1 , ..., N, ∀ t ∈ [ t n , t n +1 ] . (2.16)Now we can write the semi-discrete scheme (2.8) as ∂∂t ˜ U j ( t ) = ∂∂t ( A j ( t ) · U j ( t )) = L ( U j − ( t ) , U j ( t ) , U j +1 ( t ) , t ) , (2.17)where the spatial discretization operator on the RHS of (2.8) is denoted as L ( U j − ( t ) , U j ( t ) , U j +1 ( t ) , t ) . (4) RK time discretization and fully discrete scheme. Next, we describe the fully discretescheme with method-of-lines RK discretization of the time derivative. There are two main stepsinvolved here. 7.
Obtain the initial condition of (2.17) by an L projection of u ( x, t n ) from backgroundcells onto upstream cells I ∗ j . That is, A j ( t n ) U nj = Ç (cid:90) ˜ I j ( t n ) u ( x, t n ) ψ j, ( x, t n ) dx, · · · , (cid:90) ˜ I j ( t n ) u ( x, t n ) ψ j,k ( x, t n ) dx å T (2.18)= (cid:32)(cid:90) I (cid:63)j u ( x, t n )Ψ j, ( x + α j ∆ t n ) dx, · · · , (cid:90) I (cid:63)j u ( x, t n )Ψ j, ( x + α j ∆ t n ) dx (cid:33) T . (2.19)The integrals over the upstream cells above can be evaluated in the same fashion as the SLDG scheme [1].2. Update (2.17) from U nj to U n +1 j . We apply the SSP explicit RK methods [19] as in amethod of lines approach. In particular, the time-marching algorithm using an s -stage RKmethod follows the procedure below:(a) Get the mesh information of the dynamic element ˜ I ( l ) j , l = 0 , · · · , s on RK stages by˜ x j ± ( t ) = x j ± + ( t − t n +1 ) ν j ± .(b) For RK stages i = 1 , · · · , s , let t ( l ) = t n + d l ∆ t n , compute A j ( t ( i ) ) · U ( i ) j = i − (cid:88) l =0 î α il A j ( t ( l ) ) · U ( l ) j + β il ∆ t L Ä U ( l ) j − , U ( l ) j , U ( l ) j +1 , t ( l ) äó , (2.20)where α il and β il are related to RK methods, where we can update the coefficients U ( i ) j by inverting A j ( t ( i ) ) from the equation above. The coefficients for second, third andfourth order RK methods are provided in Table 2.1.Table 2.1: Parameters of some practical Runge-Kutta time discretizations. [19, 14]Order α il β il d l
12 12
13 1 1 0
34 14
23 12
01 0 0
12 12 -
13 13 23 13 heorem 2.1. (Mass conservation) Given a DG solution u h ( x, t n ) ∈ V kh and assuming the bound-ary condition is periodic, the proposed fully discrete GEL DG scheme with SSP RK time discretiza-tion of (2.17) 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 withunique flux at cell boundaries, as the mass conservation property of EL DG scheme [5]. We skipdetails for brevity. In this subsection, we discuss our treatment of inflow boundary conditions. We consider the lineartransport equation (2.1) with the initial condition and the inflow boundary condition ® u ( x,
0) = u ( x ) ,u ( x b , t ) = f ( t ) . (2.21)The proposed procedure for inflow boundary conditions follow steps below. At the outflow bound-ary, characteristics will go to the interior of domain, hence the original GEL DG algorithm couldbe directly applied. t n t n +1 x x x b x − x − x (cid:63) x (cid:63) x (cid:63) Ω Ω (a) t n t n +1 t (cid:63) − t (cid:63) − x x x b x − x − Ω ∗− Ω ∗ I b (b) Figure 2.2: Illustration on ghost cells intersecting inflow boundary.Step 1: Set up ghost cells. We first set up a ghost region which is sufficiently large, on whichdiscretizations are performed to define ghost cells. For example, in the 1D setting, for
CF L < x − , x − ], [ x − , x ] as illustrated in Figure 2.2 (a).Step 2: Obtain the DG solutions on ghost cells. We find DG solution at t n on ghost cells bytracking information along characteristics from boundary data in the semi-Lagrangian fashion. In9rder to do this, we first find the velocity field a ( x, t ) on ghost regions, which could be done by anatural extrapolation from interior of domain. In particular, we consider the following problem ® u t + ( a ( x, t ) u ) x = 0 , x at the ghost region u ( x b , t ) = f ( t ) , (2.22)where a ( x, t ) at the ghost region can be approximated by extrapolations from the interior of domain.As shown in Figure 2.2(b), there is Ω ∗ bounded by characteristic curves emanating from boundariesof ghost cells. We let the test function ψ ( x, t ) satisfies the adjoint problem with ∀ Ψ ∈ P k ( I b ), (cid:40) ψ t + a ( x, t ) ψ x = 0 ,ψ ( x, t n ) = Ψ( x ) , x ∈ [ x − , x ] . (2.23)Integrate ((2.22) · ψ + (2.23) · u ) over Ω ∗ , we have (cid:90) Ω ∗ ( uψ ) t + ( a ( x, t ) uψ ) x dxdt = 0 . (2.24)Using Green formula, we can get (cid:90) x x − u ( x, t n )Ψ( x ) dx = (cid:90) t ∗− t n a ( x b , t ) u ( x b , t ) ψ ( x b , t ) dt. (2.25)As in the SL DG [3, 2, 1], Gauss quadrature rule can be applied to evaluate the right-hand side ofthe above equation. Similar procedure can be used to obtain DG solutions on the ghost cells.Step 3: Update solution. Once the solution on ghost cells are available, we can update thesolution following the GEL RK DG procedure described previously. We extend the GEL RK DG algorithm to 2D problems via dimensional splitting [17]. Consider alinear 2D transport equation u t + ( a ( x, y, t ) u ) x + ( b ( x, y, t ) u ) y = 0 , ( x, y ) ∈ Ω , (2.26)with a proper initial condition u ( x, y,
0) = u ( x, y ) and boundary conditions. Here ( a ( x, y, t ) , b ( x, y, t ))is a velocity field. The domain Ω is partitioned into rectangular meshes with each computationalcell A ij = [ x i − , x i + ] × [ y j + , y j + ], where we use the piecewise Q k tensor-product polynomialspaces.1. We first locate ( k + 1) tensor-product Gaussian nodes on cell A ij : ( x i,p , y j,q ) , p, q = 0 , ..., k .For example, see Figure 2.3 (left) for the case of k = 3.10igure 2.3: Illustration of the 2D GEL RK DG scheme via Strang splitting. k = 3.2. Then, the equation (2.26) is split into two 1D advection problems based on the quadraturenodes in both x − and y − directions: u t + ( a ( x, y, t ) u ) x = 0 , (2.27) u t + ( b ( x, y, t ) u ) y = 0 . (2.28)Based on a 1D GEL RK DG formulation, the split equations (2.27) and (2.28) are evolvedvia Strang splitting over a time step ∆ t as follows. • Evolve 1D equation (2.27) at different y (cid:48) j,q s with different velocity for a half time-step∆ t/
2, see Figure 2.3 (middle). For each y j,q , the ( k + 1) point values are mapped to a P k polynomial per cell, then the 1D equation (2.27) is evolved by the proposed GEL RKDG scheme. Finally, we can map the evolved P k polynomial back to the ( k + 1) pointvalues to update the solution. • Evolve 1D equation (2.28) at different x (cid:48) i,p s for a full time-step ∆ t as above, see Figure2.3 (right). • Evolve 1D equation (2.27) at different y (cid:48) j,q s for another half time-step ∆ t/ Equivalence of semi-discrete GEL DG and EL DG methods
In this section, we first study the equivalence between the evolution step of the GEL DG and SLDG methods for a linear constant problem (2.1) in section 3.1; the equivalence between the GELDG and EL DG methods for a linear variable coefficient problem (2.1) is also presented in section3.2. In [5], we established that the evolution step in the EL DG scheme is the same as the ALEDG method, for which theoretical stability analysis are performed in [15].
The extra degree of freedom in the GEL DG scheme design, compared with the SL DG, is thespace-time partition and the adjoint problem for the test function. For a linear constant coefficientadvection equation, if we assume the exact space-time partition as the SL DG method, while varyingthe velocity field of the modified adjoint problem in GEL DG, we show below that the semi-discreteGEL DG scheme is equivalent to the SL DG scheme.
Theorem 3.1.
For linear constant coefficient equation (2.1) with a ( x, t ) = 1 , GEL DG schemewith the exact space-time partition ν j ± = 1 and a perturbation of velocity α j = 1 + c, c (cid:54) = 0 in (2.2) , is equivalent to SL DG scheme in semi-discrete form.Proof. In this case, the GEL DG scheme (2.6) is reduced into ddt (cid:90) ˜ I j ( t ) ( u GELDG ( x, t ) ψ GELDG ( x, t )) dx = − (cid:90) ˜ I j ( t ) cu GELDG ψ GELDGx dx. (3.1)We can rewrite the scheme as in integral form (cid:90) I j ( u GELDG ( x, t n +1 )Ψ( x )) dx = (cid:90) ˜ I j ( t n ) ( u GELDG ( x, t n ) ψ GELDG ( x, t n )) dx − c (cid:90) t n +1 t n (cid:90) ˜ I j ( τ ) u GELDG ( x, τ ) ψ GELDGx ( x, τ ) dxdτ = (cid:90) ˜ I j ( t n ) ( u n ( x )Ψ( x + (1 + c )∆ t )) dx − c (cid:90) t n +1 t n (cid:90) ˜ I j ( τ ) u GELDG ( x, τ )Ψ x ( x + (1 + c )( t n +1 − τ )) dxdτ. = RHS
GELDG , (3.2)where ˜ I j ( t n ) = [ x j − − ∆ t, x j + − ∆ t ] and ˜ I j ( t ) = [ x j − + t − t n +1 , x j + + t − t n +1 ] . We also12ewrite the SL DG scheme as in integral form (cid:90) I j ( u SLDG ( x, t n +1 )Ψ( x )) dx = (cid:90) ˜ I j ( t n ) ( u n ( x ) ψ SLDG ( x, t n )) dx = (cid:90) ˜ I j ( t n ) ( u n ( x )Ψ( x + ∆ t )) dx . = RHS
SLDG , (3.3)If we assume u GELDG ( x, τ ) = u SLDG ( x, τ ), τ ∈ [ t n , t n +1 ) then RHS
SLDG − RHS
GELDG = (cid:90) ˜ I j ( t n ) ( u n ( x )Ψ( x + ∆ t )) dx − (cid:90) ˜ I j ( t n ) ( u n ( x )( x, t n )Ψ( x + (1 + c )∆ t )) dx + c (cid:90) t n +1 t n (cid:90) ˜ I j ( t n ) u SLDG ( ξ + ( τ − t n ) , τ )Ψ x ( ξ + c ( t n +1 − τ ) + ∆ t ) dξdτ = (cid:90) ˜ I j ( t n ) u n ( x )(Ψ( x + ∆ t ) − Ψ( x + (1 + c )∆ t )) dx − (cid:90) ˜ I j ( t n ) u n ( ξ ) (cid:90) t n +1 t n − c Ψ x ( ξ + c ( t n +1 − τ ) + ∆ t ) dτ dξ = (cid:90) ˜ I j ( t n ) u n ( x )(Ψ( x + ∆ t ) − Ψ( x + (1 + c )∆ t )) dx − (cid:90) ˜ I j ( t n ) u n ( ξ ) (cid:90) t n +1 t n Ψ τ ( ξ + c ( t n +1 − τ ) + ∆ t ) dτ dξ = 0 , This verifies that the u GELDG = u SLDG for semi-discrete schemes.The fully discrete GEL RK DG and SL DG scheme are not equivalent for any P k approximationspaces. The equivalence holds true for the special case of GEL RK DG P and P schemes, asspecified in the Theorem below. Theorem 3.2.
Under the same condition as Theorem 3.1, the fully discrete GEL RK DG and SLDG scheme are equivalent for P k ( k ≤ ) approximation spaces with any RK time discretization.Thus GEL RK DG schemes with P k ( k ≤ ) approximation spaces are unconditionally stable.Proof. We first consider the GEL DG method with forward-Euler time discretization by (3.1) (cid:90) I j u GELDG ( x, t n +1 )Ψ( x ) dx = (cid:90) ˜ I j ( t n ) u n ( x )Ψ( x + (1 + c )∆ t ) dx − c ∆ t (cid:90) ˜ I j ( t n ) u n ( x )Ψ x ( x + (1 + c )∆ t ) dx. = RHS GELDG . We compute the right hand side
RHS
SLDG − RHS GELDG = (cid:90) ˜ I j ( t n ) u n ( x )(Ψ( x + ∆ t ) − Ψ( x + (1 + c )∆ t )) dx + ∆ t (cid:90) ˜ I j ( t n ) cu n ( x )Ψ x ( x + (1 + c )∆ t ) dx = (cid:90) ˜ I j ( t n ) u n ( x )(Ψ( x + ∆ t ) − Ψ( x + (1 + c )∆ t ) + c ∆ t Ψ x ( x + (1 + c )∆ t )) dx x ) ∈ P k , k ≤
1, then
RHS
SLDG − RHS GELDG = 0. That is, the GEL DG and SL DGscheme are equivalent for P k , k ≤ Remark 3.3.
It can be shown that under the same assumption as Theorem 3.1 with P approxi-mation space, the fully discrete GEL RK DG method, when the SSP RK2 and RK3 are applied fortime discretization in Table 2.1, is equivalent to the SL DG. Thus the scheme is unconditionallystable. However, the GEL RK DG, when coupled with forward Euler time discretization, is notequivalent to the SL DG and is not unconditionally stable. Theorem 3.4.
For a linear transport problem with variable coefficient (2.1) , GEL DG scheme isequivalent to EL DG scheme in the semi-discrete form, assuming that they have the same space-timepartition.Proof.
We first consider EL DG scheme for scaler equation (2.1), which is formulated as ddt (cid:90) ˜ I j ( t ) ( uψ ) dx = − Ä ˆ F ψ ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + Ä ˆ F ψ ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t ) + (cid:90) ˜ I j ( t ) F ψ x dx, (3.4)with F ( u ) . = ( a − α ) u , α ( x, t ) = − ν j − x − ˜ x j + ( t )∆ x j ( t ) + ν j + x − ˜ x j − ( t )∆ x j ( t ) ∈ P ([˜ x j − ( t ) , ˜ x j + ( t )]) , (3.5)and ˆ F is a monotone numerical flux. We represent u in the form of (2.14), and take ψ as our basisfunction ˜ ψ j,m ( x, t ) satisfy (2.13) with Ψ j,m ( x ) being orthonormal basis of the space P k ( I j ), whichfollows the adjoint problem of EL DG: ® ψ t + α ( x, t ) ψ x = 0 , ( x, t ) ∈ Ω j ,ψ ( t = t n +1 ) = Ψ( x ) , ∀ Ψ( x ) ∈ P k ( I j ) . (3.6)Then we have (cid:90) ˜ I j ( t ) ˜ ψ j,i ( x, t ) ˜ ψ j,m ( x, t ) dx = ∆ x j ( t ) (cid:90) I j Ψ j,i ( x )Ψ j,m ( x ) dx = δ i,m ∆ x j ( t ) . We can write the EL DG scheme (3.4) into a matrix form ddt ( U j ( t )∆ x j ( t )) = − Ä ˆ F ˜ ψ j ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + Ä ˆ F ˜ ψ j ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t ) + ˜ B j ( t ) U j ( t ) , (3.7)14here ˜ B j ( t ) is a time-dependent matrix consisting of [ ˜ B j ( t )] i,m = (cid:82) ˜ I j ( t ) ( a ( x, t ) − α ) ˜ ψ j,m ( ˜ ψ j,i ) x dx, i =0 , ...k, m = 0 , ...k , and ˜ ψ j = ( ˜ ψ j, , ..., ˜ ψ j,k ) T is a vector of size ( k + 1) × ddt (cid:90) ˜ I j ( t ) ( uψ j,i ( x, t )) dx = − Ä ˆ F ψ j,i ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + Ä ˆ F ψ j,i ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t ) + (cid:90) ˜ I j ( t ) ( a ( x, t ) − α ) u ( ψ j,i ) x dx, (3.8)where ψ j,i is the test function as defined in (2.11). Then, we rewrite the GEL DG scheme (3.8)into the matrix form ddt ( A j ( t )∆ x j ( t ) U j ( t )) = − Ä ˆ F ψ j ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + Ä ˆ F ψ j ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t ) + B j ( t ) U j ( t ) , (3.9)where [ A j ( t )] i,m = (cid:82) ˜ I j ( t ) ˜ ψ j,m ( x, t ) ψ j,i ( x, t ) dx/ ∆ x j ( t ), [ B j ( t )] i,m = (cid:82) ˜ I j ( t ) ( a ( x, t ) − α j ) ˜ ψ j,m ( ψ j,i ) x dx and ψ j = ( ψ j, , ..., ψ j,k ) T . Applying the product rule to the LHS of (3.9), and with manipulationsof the equation, we have ddt (∆ x j ( t ) U j ( t )) = A − j ( t ) (cid:16) − Ä ˆ F ψ j ä (cid:12)(cid:12)(cid:12) j + + Ä ˆ F ψ j ä (cid:12)(cid:12)(cid:12) j + + B j ( t ) U j ( t ) − ˙ A j ( t )∆ x j ( t ) U j ( t ) (cid:17) . (3.10)As { ˜ ψ j,i ( x, t ) } ki =0 and { ψ j,i ( x, t ) } ki =0 are basis of P k ( ˜ I j ( t )), we can represent ψ j,i ( x, t ) as ψ j,i ( x, t ) = k (cid:88) n =0 ( ψ j,i ( x, t ) , ˜ ψ j,n ( x, t )) ˜ I j ( t ) ˜ ψ j,n ( x, t ) / ∆ x j ( t ) = [ A j ( t )] i ˜ ψ j , where [ A j ( t )] i is the i th row of matrix A j ( t ). That is, A − j ( t ) ψ j = ˜ ψ j , (3.11)which leads to A − j ( t )( − Ä ˆ F ψ j ä (cid:12)(cid:12)(cid:12) j + + Ä ˆ F ψ j ä (cid:12)(cid:12)(cid:12) j + ) = − Ä ˆ F ˜ ψ j ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t ) + Ä ˆ F ˜ ψ j ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t ) . Next, we compute ˙ A j ( t ), by ∆ x j ( t ) = ∆ x +( ν j − − ν j + )( t n +1 − t ) and ∂ x ( α ) = ( ν j + − ν j − ) / ∆ x j ( t )[ ˙ A j ( t )] i,m = 1∆ x j ( t ) ddt (cid:90) ˜ I j ( t ) ˜ ψ j,m ( x, t ) ψ j,i ( x, t ) dx + ( ν j − − ν j + )(∆ x j ( t )) (cid:90) ˜ I j ( t ) ˜ ψ j,m ( x, t ) ψ j,i ( x, t ) dx = 1∆ x j ( t ) (cid:90) ˜ I j ( t ) Ä ∂ t ( ˜ ψ j,m ) ψ j,i + ˜ ψ j,m ∂ t ( ψ j,i ) + ∂ x ( α ˜ ψ j,m ψ j,i ) ä dx + ( ν j − − ν j + )(∆ x j ( t )) (cid:90) ˜ I j ( t ) ˜ ψ j,m ( x, t ) ψ j,i ( x, t ) dx = 1∆ x j ( t ) (cid:90) ˜ I j ( t ) Ä − α∂ x ( ˜ ψ j,m ) ψ j,i − ˜ ψ j,m α j ∂ x ( ψ j,i ) + ∂ x ( α ) ˜ ψ j,m ψ j,i + α∂ x ( ˜ ψ j,m ) ψ j,i + α ˜ ψ j,m ∂ x ( ψ j,i ) ä dx + ( ν j − − ν j + )(∆ x j ( t )) (cid:90) ˜ I j ( t ) ˜ ψ j,m ( x, t ) ψ j,i ( x, t ) dx
15 1∆ x j ( t ) Ç (cid:90) ˜ I j ( t ) ( α − α j ) ˜ ψ j,m ∂ x ( ψ j,i ) dx + ( ν j + − ν j − )∆ x j ( t ) (cid:90) ˜ I j ( t ) ˜ ψ j,m ψ j,i dx å + ( ν j − − ν j + )(∆ x j ( t )) (cid:90) ˜ I j ( t ) ˜ ψ j,m ( x, t ) ψ j,i ( x, t ) dx = 1∆ x j ( t ) (cid:90) ˜ I j ( t ) ( α − α j ) ˜ ψ j,m ( x, t ) ∂ x ( ψ j,i )( x, t ) dx. So we have [ B j ( t )] i,m − [ ˙ A j ( t )] i,m ∆ x j ( t ) = (cid:82) ˜ I j ( t ) ( a ( x, t ) − α ) ˜ ψ j,m ( x, t ) ∂ x ( ψ j,i )( x, t ) dx . Then we caneasily get A − j ( t )( B j ( t ) − ˙ A j ( t )∆ x j ( t )) = ˜ B j ( t ) by (3.11), which shows the equivalence of (3.7) and(3.10). Remark 3.5.
For general nonlinear problems, the equivalence of semi-discrete GEL DG and ELDG methods can be established in a similar way, given the same space-time partition.
Remark 3.6. (Fully discrete case) The fully discrete EL RK DG scheme is known as equivalent tothe fully discrete ALE DG scheme combined with one extra step of solution projection (which doesnot affect stability). Thus the stability result of ALE DG method for linear conservation laws [28]can be directly applied to assess the stability property of fully discrete EL RK DG scheme. Thestability of fully discrete GEL RK DG method is still theoretically open, and is being investigatednumerically in the numerical section.
We consider the property of DGCL, which requires that the numerical scheme reproduces exactlya constant solution under the geometric parameters of numerical schemes. As shown in [8], sat-isfaction of DGCL is a necessary and sufficient condition for a numerical scheme to preserve thenonlinear stability of its fixed grid counterpart. In the context of GEL DG scheme, the geomet-ric parameters refer to the parameters involved in the space-time partition on which the PDE isevolved. It was established in [15] that ALE-DG scheme (the evolution step in the EL DG method)satisfies the DGCL for 1D problems, and for high D problems with the time integrator which holdsthe accuracy not less than the value of the spatial dimension [10]. Below we shown the conditionsunder which the DGCL holds for GEL DG method when coupled with forward Euler time dis-cretization, see Table 4.2. As a direct consequence, we find the DGCL no longer holds for generalGEL RK DG schemes. Proposition 4.1 is numerically verified in Table 4.3.16able 4.2: The related parameters settings about RHS and DGCL of GEL DG method for u t + au x =0 . Ψ RHS condition for DGCL1 ∆ x – x ∆ t ∆ x ( ν j − − ν j + ) + ∆ t ∆ x ( ν j − − ν j + )( ν j − + ν j + 12 − α j ) ν j + = ν j − x − − ∆ t x [( ν j + − α j ) + ( ν j − − α j ) ] + t x [( ν j + − α j ) − ( ν j − − α j ) ] ν j + = ν j − = α j Table 4.3: DGCL. The error (cid:107) u h ( x, t ) − u ( x, t ) (cid:107) L ∞ of the GEL DG scheme solving u t + u x = 0 withinitial condition u ( x,
0) = 1, T = 1, CF L = 0 .
1. We test schemes with different choices of ν j ± asparameters for space-time partition and α j as parameters for adjoint problems. ν j ± α j . x sin( x j ) P − − − P − − − P − − − P − − − P − − − P − − − P − − − x sin( x j ± ) P − − − P − − − Proposition 4.1.
Under the conditions specified in Table 4.2 that the GEL DG method coupledwith forward Euler time discretization satisfies the DGCL.
Proof.
For the GEL DG scheme with Forward-Euler time discretization, we can get (cid:90) ˜ I j ( u n +1 j ψ ( x, t n +1 )) dx = (cid:90) I ∗ j ( u nj ψ ( x, t n )) dx + ∆ t ( − Ä ˆ F ψ ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j + 12 ( t n ) + Ä ˆ F ψ ä (cid:12)(cid:12)(cid:12)(cid:12) ˜ x j − ( t n ) )+ ∆ t (cid:90) I ∗ j ( a − α j ) uψ x dx . = RHS. (4.1)We let u nj = 1, RHS = (cid:90) I ∗ j Ψ( x + α j ∆ t ) dx + ∆ t ( − ( a − ν )Ψ( x + α j ∆ t ) | ˜ x j + 12 ( t n ) + ( a − ν )Ψ( x + α j ∆ t ) | ˜ x j − ( t n ) + ( a − α j )Ψ( x + α j ∆ t ) | ˜ x j + 12 ( t n ) − ( a − α j )Ψ( x + α j ∆ t ) | ˜ x j − ( t n ) )= (cid:90) I ∗ j Ψ( x + α j ∆ t ) dx + ∆ t ( − ( α j − ν )Ψ( x + α j ∆ t ) | ˜ x j + 12 ( t n ) + ( α j − ν )Ψ( x + α j ∆ t ) | ˜ x j − ( t n ) ) . With different choices of Ψ, we show conditions in Table 4.2, under which the DGCL is satisfied.From the above proposition, the GEL DG method with P polynomial space and forward Eulersatisfies DGCL. However, for the method with high order polynomial spaces, DGCL fails for generalGEL RK DG schemes. On the other hand, MPP is the principle that numerical solutions will be17ounded by the maximum and minimum of initial condition. One can show that if a scheme satisfiesthe MPP, it automatically satisfies the DGCL. Below, we are going to first show that first orderGEL DG scheme satisfies the MPP property. Then we apply two MPP limiters to high order GELRK DG schemes to get the MPP, hence the DGCL property. Proposition 4.2. (MPP property of the first order GEL DG scheme.) Let u m = min u ( x ) , u M =max u ( x ), then the first order GEL DG solution ¯ u nj ∈ [ u m , u M ], ∀ j, n under the condition ∆ t ≤ ∆ xα where α = max j | a − ν | j + . Proof.
The first order GEL DG scheme reads (cid:90) ˜ I j ¯ u n +1 j dx = (cid:90) I ∗ j ¯˜ u nj dx − ∆ t Ä ˆ h j + − ˆ h j − ä where ˆ h is the first order monotone flux, e.g. the Lax-Friedrich flux,ˆ h j + = ˆ F | ˜ x j + 12 ( t n ) = ( a − ν ) j + (¯˜ u nj +1 + ¯˜ u nj )2 − α (¯˜ u nj +1 − ¯˜ u nj )2 , and ˜ u nj is the reaveraging of solution u nj on I ∗ j . ¯˜ u nj ∈ [ u m , u M ], since u nj ∈ [ u m , u M ]. Let λ = ∆ t ∆ x ,we get ¯ u n +1 j = ∆ x ∗ j ∆ x ¯˜ u nj − λ [ˆ h j + − ˆ h j − ]= (1 − λα )¯˜ u nj + λ α − ( a − ν ) j + u nj +1 + λ ( a − ν ) j − + α u nj − . Under the condition of ∆ t ≤ ∆ xα with α = max j | a − ν | j + , the coefficients 1 − λα , α ± ( a − ν ) j − are all positive, hence ¯ u n +1 j is a convex combination of ¯˜ u nj , ¯˜ u nj +1 , ¯˜ u nj − . Therefore, ¯ u n +1 j ∈ [ u m , u M ].Next we propose to apply MPP limiters to high order GEL RK DG schemes to obtain the MPPproperty, leading to DGCL in a general setting. There are two MPP limiters that we will make useof: one is the rescaling limiter for DG polynomials [26, 27], and the other is the parametrized fluxlimiter in a sequence of work [22, 24, 23]. While the rescaling limiter works very well for RKDGschemes, it often leads to extra time step restrictions for MPP property and reduction in temporalorder of convergence [26]. For the GEL RK DG scheme, we propose to apply the rescaling limiter[26, 27] only for solutions at t n and in the solution remapping step, and apply the flux limiter only18n the final stage of RK methods [22]. Neither the rescaling limiter nor the flux limiter are appliedin internal stages of RK methods to avoid the temporal order reduction of RK methods.Below, we demonstrate the procedure of the proposed MPP flux limiter in the context of a thirdorder SSP RK time discretization.Step 1: We first use the polynomial rescaling MPP limiter [26] to preserve the MPP property forthe DG polynomial u ( x, t n ) on I j , then perform the L projection of u ( x, t n ) from background cellsonto upstream cells I ∗ j . We assume ˜ u nj as the reprojected polynomial on I ∗ j , then ¯˜ u nj ∈ [ u m , u M ].Step 2: We propose to apply the parametrized MPP flux limiters [23] in the context of a movingmesh to guarantee ¯ u n +1 j ∈ [ u m , u M ]. Taking ψ = 1 in (2.8), we have ddt (cid:90) ˜ I j ( t ) u j ( x, t ) dx = ˆ H j + − ˆ H j − , (4.2)where ˆ H is the high order monotone flux. With the third order SSP RK method, the update ofcell averages in equation (4.2) can be written as¯ u n +1 j = ∆ x ∗ j ∆ x ¯˜ u nj − λ [ ˆ H rkj + − ˆ H rkj − ] , where ˆ H rkj + = ˆ H nj + + ˆ H (1) j + + ˆ H (2) j + with ˆ H ( i ) , i = 0 , , u h ( x, t ( i ) ) at each RK stage. The MPP flux limiter is proposed to replace the numerical flux ˆ H rkj + by a modified one ˜ H rkj + = θ j + ( ˆ H rkj + − ˆ h j + ) + ˆ h j + , where ˆ h j + is the numerical flux from afirst order scheme as in Proposition 4.2 and θ j + ∈ [0 , ∀ j . The parameter θ j + is set to be asclose to 1 as possible, and to ensure ¯ u n +1 j ∈ [ u m , u M ], for which sufficient inequalities to have are, ∀ j , λθ j − ( ˆ H rkj − − ˆ h j − ) − λθ j + ( ˆ H rkj + − ˆ h j + ) − Γ Mj ≤ , (4.3) λθ j − ( ˆ H rkj − − ˆ h j − ) − λθ j + ( ˆ H rkj + − ˆ h j + ) − Γ mj ≥ , (4.4)with Γ Mj = u M − ∆ x ∗ j ∆ x ¯ u nj + λ (ˆ h j + − ˆ h j − ) , Γ mj = u m − ∆ x ∗ j ∆ x ¯ u nj + λ (ˆ h j + − ˆ h j − ) . Here we adjust the terms Γ Mj and Γ mj from [23] in the context of moving meshes, and obtain theparameter θ j + satisfying (4.3) and (4.4) for all j, thus guarantee u n +1 h ∈ [ u m , u M ]. Note thatsuch θ always exists, since θ j + = 0 is a solution to (4.3) and (4.4) for all j . Then we go back toStep 1 to apply the polynomial rescaling MPP limiter again to ensure DG polynomials u h ( x, t n +1 )satisfy the MPP property. 19e call the GEL RK DG method with the above described limiter ‘GELDGMPPlimiter’. Wealso apply the polynomial rescaling MPP limiter [26] at each RK internal stages, for which thescheme is termed ‘zhangMPPlimiter’. We will compare numerical performance, in terms of errorand time stepping size for numerical stability, of these two limiters in the next section. Remark 4.3.
For a linear variable coefficient equation, MPP property is lost, but the PP propertystays valid. A PP limiter can be applied in a similar fashion to preserve the PP property.
In this section, we perform numerical experiments for linear transport problems, where we set thetime stepping size as ∆ t = CF La ∆ x for 1D and ∆ t = CF L a ∆ x + b ∆ y , where a and b are maximum transportspeed in x - and y -directions respectively. We mainly study the following aspects: the spatial order ofconvergence by using small enough time stepping size, the temporal order of convergence by varyingCFL, numerical stability under a large time stepping size. When applicable, we also present theEL DG solutions [5] for comparison. Example 5.1. (1D linear transport equation with constant coefficient.) We consider a simple 1Dtransport equation u t + u x = 0 , x ∈ [0 , π ] , (5.1)with the smooth initial data u ( x,
0) = sin( x ) and exact solution u ( x, t ) = sin( x − t ). For the constantcoefficient problem, the proposed GEL DG method, if using the exact velocity field for space-timepartition and the adjoint problem, is the same as EL DG and SL DG. Here we perturb the velocityat cell boundaries, i.e. ν j + and the velocity in the modified adjoint problem i.e. α j in (2.2) to getGELDG1, GELDG2 and GELDG3 schemes respectively. Parameters of these GEL DG methods aregiven in Table 5.4. Table 5.5 reports the spatial accuracies of these methods for this example withTable 5.4: The numerical parameters of ELDG, GELDG1, GELDG2 and GELDG3 method for(5.1). ELDG GELDG1 GELDG2 GELDG3mesh ν j + x j + )∆ x x j + )∆ x x j + )∆ x adjoint α j – 1 1 + sin( x j )∆ x x j )∆ x the same time stepping size. The proposed GEL DG methods are found to perform comparably20s the ELDG method. We vary time stepping size, with fixed well-resolved spatial meshes, andplot error vs. CF L in Figure 5.4 for these schemes with P -SSP RK2 solutions (left) and P -SSPRK3 solutions (right) at a long integration time T = 100. GELDG2 methods are found to beunconditionally stable with the space-time partition exactly following the characteristics, which isconsistent with Theorem 3.2 and Remark 3.3. GELDG3 behaves closer to ELDG and performsbetter than GELDG1, as the adjoint problem and the space-time partition are more closely related.It indicates that, designing the GEL DG scheme associating the adjoint problem with the space-timepartition is advantageous for better performance of the scheme. In particular, the adjoint problemis determined by the space-time partition in the EL DG algorithm for which best numerical stabilityis observed and theoretically investigated in [5].Further, we apply higher order RK methods for time discretization, as we are interested inusing relatively large time stepping size. We show error vs. CF L in Figure 5.5 for GELDG1and GELDG3 schemes P with SSP RK2, RK4 (left) and P SSP RK3, RK4 (right) at a longintegration time T = 100. We can concluded that GELDG3 is better than GELDG1 in terms ofstability and the higher order RK help with reducing the error magnitute when large time steppingsize is used. We note that in both Figure 5.4 and 5.5, the CFL allowed with stability is much largerthan that of the RK DG method which is k +1 .Table 5.5: 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 . CF L = 0 .
18 for all P and P schemes,respectively. Mesh L error Order L error Order L error Order L error Order P ELDG P GELDG1 P GELDG2 P GELDG340 6.08E-04 – 6.08E-04 – 6.37E-04 – 6.08E-04 –80 1.55E-04 1.97 1.55E-04 1.97 1.59E-04 2.00 1.55E-04 1.97160 3.84E-05 2.02 3.84E-05 2.02 3.90E-05 2.03 3.84E-05 2.02320 9.77E-06 1.98 9.77E-06 1.98 9.83E-06 1.99 9.77E-06 1.98 P ELDG P GELDG1 P GELDG2 P GELDG340 7.69E-06 – 2.60E-05 – 7.25E-06 – 7.71E-06 –80 9.45E-07 3.03 1.91E-06 3.77 9.23E-07 2.97 9.45E-07 3.03160 1.18E-07 3.00 1.67E-07 3.51 1.17E-07 2.98 1.18E-07 3.00320 1.41E-08 3.07 1.63E-08 3.36 1.40E-08 3.06 1.41E-08 3.07
Finally, we verify the DGCL property of the scheme when the proposed MPP limiters areapplied. Table 5.6 shows that we regain the DGCL property by applying ‘GELDGMPPlimiter’to GEL RK DG schemes, whose DGCL property are assessed in Table 4.3. Next, we present inFigure 5.6 the L ∞ error versus CF L of GEL DG methods with ‘GELDG MPP limiter’ and ‘zhang21igure 5.4: The L ∞ error versus CF L of ELDG methods, GELDG1, GELDG2 and GELDG3methods P (left) with SSP RK2 and P (right) with SSP RK3 time discretization for (5.1) withinitial condition u ( x,
0) = sin( x ). A long time simulation is performed with T = 100 and mesh size N = 160.Figure 5.5: The L ∞ error versus CF L of GELDG1 and GELDG3 methods P (left) with SSP RK2,RK4 and P (right) schemes SSP RK3, RK4 time discretization for (5.1) with initial condition u ( x,
0) = sin( x ). A long time simulation is performed with T = 100 and mesh size N = 160.MPP limiter’ P and P schemes. We apply the RK4 for time integration. From Figure 5.6, wecan observe better stability and accuracy for GEL DG method with ‘GELDGMPPlimiter’ which isapplied only in the final RK stage, compared with one with ‘zhangMPPlimiter’ which is applied inevery RK intermediate stages, especially for schemes with P polynomial space. We also test theMPP property for a step function initial condition: u ( x,
0) = , < x < , , otherwise. (5.2)The computational domain is [0,90]. The solution of GELDG2 method for P with CF L = 1, N = 160 and with RK4 time discretization without (left) and with (right) ‘GELDGMPPlimiter’are shown in Figure 5.7. We can observe the MPP property for GEL DG with ‘GELDGMPPlimiter’.22able 5.6: Compute u h ( x, t ) − u ( x, t ) to check if all P , P and P GEL DG schemes with‘GELDGMPPlimiter’ satisfy DGCL for (5.1) with initial condition u ( x,
0) = 1, T = 1, CF L = 0 . ν j ± α j . x sin( x j ) P − − − P − − − P − − − P − − − P − − − P − − − P − − − x sin( x j ± ) P − − − P − − − Figure 5.6: The L ∞ error versus CF L of GELDG1, GELDG2 and GELDG3 methods with ‘GELDGMPP limiter’ and ‘zhang MPP limiter’ for P (left) and P (right) with RK4 time discretizationfor (5.1) with initial condition u ( x,
0) = sin( x ). A long time simulation is performed with T = 100and mesh size N = 160.Figure 5.7: The solution u of GELDG2 method P without (left) and with (right) ‘GELDGMP-Plimiter’ with RK4 time discretization and CF L = 1 for (5.1) with initial condition (5.2). A longtime simulation is performed with T = 40 and mesh size N = 160.23 xample 5.2. (1D transport equation with variable coefficients.) Consider u t + (sin( x ) u ) x = 0 , x ∈ [0 , π ] (5.3)with initial condition u ( x,
0) = 1 and the periodic boundary condition. The exact solution is givenby u ( x, t ) = sin(2 tan − ( e − t tan( x )))sin( x ) . (5.4)The related parameters settings are given in Table 5.7. The expected spatial convergence of ELDGand GELDG are shown in Table 5.8. In Figure 5.8, we plot the L ∞ error versus CF L of ELDGand GELDG scheme with P (left) and P (right) polynomial spaces. The following observationsare made: (1) both methods perform similarly around and before CF L = 1, which is well abovethe stability constraint of the RK DG method 1 / (2 k + 1); (2) after CF L = 1 and before stabilityconstraint of the method, the temporal convergence order is observed to be consistent with theorder of RK discretization; (3) EL RK DG has better performance than GEL RK DG with thesame mesh.In addition, we test the proposed inflow boundary condition for the following problem u t + (sin( x ) u ) x = 0 , x ∈ [ π , π ] u ( x,
0) = 1 ,u ( π , t ) = sin(2 tan − ( e − t tan( π ))) . (5.5)We use the GEL DG method with PP limiter to solve this problem, the L and L ∞ errors for P and P are shown in Table 5.9. The optimal convergence rate are observed. Besides, we showthe L ∞ error versus CF L of GEL DG method with PP limiter in Figure 5.9. Expected temporalconvergence is observed.Table 5.7: The related parameters settings of ELDG, GELDG method for u t + (sin( x ) u ) x = 0 . setting ELDG GELDGmesh ν j + sin( x j + ) sin( x j + )adjoint α j – sin( x j ) Example 5.3. (Swirling deformation flow). Consider u t − (cos ( x y ) g ( t ) u ) x + (sin( x ) cos ( y g ( t ) u ) y = 0 , ( x, y ) ∈ [ − π, π ] (5.6)24able 5.8: 1D transport equation with variable coefficients. u t + (sin( x ) u ) x = 0 with the initialcondition u ( x,
0) = 1. T = 1. We use CF L = 0 . CF L = 0 .
18 for all P (RK2) and P (RK3)schemes, respectively. Mesh L error Order L error Order P GELDG P GELDG40 1.36E-03 – 1.36E-03 –80 3.57E-04 1.93 3.57E-04 1.93160 8.95E-05 1.99 8.95E-05 1.99320 2.31E-05 1.95 2.31E-05 1.95 P ELDG P GELDG40 5.15E-05 – 5.20E-05 –80 6.33E-06 3.03 6.37E-06 3.03160 7.84E-07 3.01 7.89E-07 3.01320 9.60E-08 3.03 9.71E-08 3.02
Figure 5.8: The L ∞ error versus CF L of ELDG methods and GELDG methods for (5.3) with theinitial condition u ( x,
0) = 1. T = 1. ∆ t = CF L ∆ x .Table 5.9: 1D transport equation with variable coefficients. u t + (sin( x ) u ) x = 0 with the initialcondition u ( x,
0) = 1 and inflow boundary condition. T = 1. We use CF L = 0 . CF L = 0 . P (RK2) and P (RK3) schemes, respectively. Mesh L error Order L error Order P GELDGwithout PP limiter with PP limiter20 5.06E-03 – 4.97E-03 –40 1.36E-03 1.90 1.36E-03 1.8880 3.55E-04 1.94 3.55E-04 1.94160 8.90E-05 1.99 8.90E-05 1.99320 2.30E-05 1.95 2.30E-05 1.95 P GELDGwithout PP limiter with PP limiter20 4.16E-04 – 4.16E-04 –40 5.20E-05 3.00 5.20E-05 3.0080 6.38E-06 3.03 6.38E-06 3.03160 7.90E-07 3.01 7.90E-07 3.01320 9.72E-08 3.02 9.72E-08 3.02 where g ( t ) = cos( πtT ) π and T = 1 .
5. The initial condition is the following smooth cosine bells (with C u ( x, y,
0) = r b cos ( r b r b π ) , if r b < r b , , otherwise, (5.7)25igure 5.9: The L ∞ error versus CF L of GELDG method with PP limiter for (5.3) with the initialcondition u ( x,
0) = 1 and inflow boundary. T = 1. ∆ t = CF L ∆ x .where r b = 0 . π , and r b = » ( x − x b ) + ( y − y b ) denotes the distance between ( x, y ) and thecenter of the cosine bell ( x b , y b ) = (0 . π, t = T /
2, then goes back to its initial shape at t = T as theflow reverses. If this problem is solved up to T , we call such a procedure one full evolution. Wetest accuracy for Q k EL DG, GEL DG methods without and with PP limiter with 4th RK and4th splitting method for k = 1 , CF L = 2 . k + 1)th order convergence is observed for these methods. Weplot the L ∞ error versus CF L of EL RK DG, GEL DG methods with Q (left) and Q (right)polynomial spaces for this case in Figure 5.10 which show that these three methods have similarstability for higher order discretization. We also plot the L ∞ error versus CF L of GEL DG methodswith PP limiter in Figure 5.11, in comparison to the one without limiter in the previous figure.Figure 5.10: The L ∞ error versus CF L of EL DG methods and GEL DG methods for (5.6) withthe smooth cosine bells (5.7) at T = 1.5. 26able 5.10: Swirling deformation flow. Q k EL DG and GEL DG methods without and withPPlimiter (k = 1, 2) for (5.6) with the smooth cosine bells (5.7) at T = 1.5.
CF L = 2 . N Q Q L -error order L ∞ -error order L -error order L ∞ -error orderEL DG 20 Figure 5.11: The L ∞ error versus CF L of GEL DG with PP limiter methods for (5.6) with thesmooth cosine bells (5.7) at T = 1.5.Lastly, we test two schemes on the swirling deformation flow (5.6) with the following setting:the computational domain is [ − π, π ] with the periodic boundary conditions and an initial conditionplotted in Figure 5.12, which consists of a slotted disk, a cone as well as a smooth hump, similar tothe one used in [16]. It is a challenging test for controlling oscillations around discontinuities. Weadopt a simple TVB limiter with M = 15 in [7] for all schemes. We simulate this problem after onefull revolutions and report the numerical solutions in Figure 5.13. For better comparison, we plot1D cuts and zoom in of the numerical solutions along with the exact solution to demonstrate theeffectiveness of the PP limiter in Figure 5.14. It is found that oscillations are well controlled withthe TVB limiter and are positivity preserving with the PP limiter. Solutions with larger CFL areobserved to dissipate less than solutions from smaller CFL.27igure 5.12: Plots of the initial profile. The mesh of 400 ×
400 is used.Figure 5.13: Plots of the numerical solutions of GEL DG schemes with TVB and PP limitersfor solving (5.6) with initial data plotted in Fig. 5.12. The final integration time T is 1.5. Themesh of 100 ×
100 is used. Left: Q GELDG-split+TVB+PPlimiter with
CF L = 2 .
2. Right: Q GELDG-split+TVB+PPlimiter with
CF L = 5 . In this paper, we develop a generalized Eulerian-Lagrangian (GEL) discontinuous Galerkin (DG)method for linear transport problems. Inflow boundary treatment is discussed. The method has theadvantages in stability under large time stepping sizes, and in mass conservation, compactness and28igure 5.14: Plots of 1D cuts and zoom in of the numerical solution for GEL DG methods withand without PP limiter for solving (5.6) with initial data Fig. 5.12. The mesh of 100 ×
100 is used.Left: numerical solution at x = 0 + π/ y = π/ π/ References [1] X. Cai, W. Guo, and J.-M. Qiu. A high order conservative semi-Lagrangian discontinuousGalerkin method for two-dimensional transport simulations.
Journal of Scientific Computing ,73(2-3):514–542, 2017.[2] X. Cai, W. Guo, and J.-M. Qiu. A high order semi-Lagrangian discontinuous Galerkin method29or Vlasov-Poisson simulations without operator splitting.
Journal of Computational Physics ,354:529–551, 2018.[3] X. Cai, W. Guo, and J.-M. Qiu. Comparison of semi-Lagrangian discontinuous Galerkinschemes for linear and nonlinear transport simulations.
Communications on Applied Mathe-matics and Computation , pages 1–31, 2020.[4] X. Cai, J. Qiu, and J.-M. Qiu. A conservative semi-Lagrangian HWENO method for theVlasov equation.
Journal of Computational Physics , 323:95–114, 2016.[5] X. Cai, J.-M. Qiu, and Y. Yang. An Eulerian-Lagrangian discontinuous Galerkinmethod for transport problems and its application to nonlinear dynamics. arXiv preprintarXiv:2002.02930 , 2020.[6] M. Celia, T. Russell, I. Herrera, and R. Ewing. An Eulerian-Lagrangian localized adjointmethod for the advection-diffusion equation.
Advances in water resources , 13(4):187–206,1990.[7] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finiteelement method for conservation laws II: general framework.
Mathematics of Computation ,52(186):411–435, 1989.[8] C. Farhat, P. Geuzaine, and C. Grandmont. The discrete geometric conservation law and thenonlinear stability of ALE schemes for the solution of flow problems on moving grids.
Journalof Computational Physics , 174(2):669–694, 2001.[9] E. Forest and R. D. Ruth. Fourth-order symplectic integration.
Physica D: Nonlinear Phe-nomena , 43(1):105–117, 1990.[10] P. Fu, G. Schn¨ucke, and Y. Xia. Arbitrary Lagrangian-Eulerian discontinuous Galerkin methodfor conservation laws on moving simplex meshes.
Mathematics of Computation , 88:2221–2255,2019.[11] E. Hairer, C. Lubich, and G. Wanner.
Geometric numerical integration , volume 31 of
Springer Series in Computational Mathematics . Springer-Verlag, Berlin, second edition, 2006.Structure-preserving algorithms for ordinary differential equations.3012] X. Hong and Y. Xia. Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for hyper-bolic equations involving δ -singularities. SIAM Journal on Numerical Analysis , 58(1):125–152,2020.[13] C.-S. Huang, T. Arbogast, and J. Qiu. An Eulerian–Lagrangian WENO finite volume schemefor advection problems.
Journal of Computational Physics , 231(11):4028–4052, 2012.[14] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted ENO schemes.
Journal ofComputational Physics , 126(1):202–228, 1996.[15] C. Klingenberg, G. Schn¨ucke, and Y. Xia. Arbitrary Lagrangian-Eulerian discontinuousGalerkin method for conservation laws: analysis and application in one dimension.
Mathe-matics of Computation , 86(305):1203–1232, 2017.[16] R. LeVeque. High-resolution conservative algorithms for advection in incompressible flow.
SIAM Journal on Numerical Analysis , 33(2):627–665, 1996.[17] 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(23):8386–8409, 2011.[18] T. Russell and M. Celia. An overview of research on Eulerian-Lagrangian localized adjointmethods (ELLAM).
Advances in water resources , 25(8):1215–1231, 2002.[19] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturingschemes.
Journal of Computational Physics , 77(2):439–471, 1988.[20] 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(1):120–163, 1999.[21] K. Wang, H. Wang, and M. Al-Lawatia. An Eulerian-Lagrangian discontinuous Galerkinmethod for transient advection-diffusion equations.
Numerical Methods for Partial DifferentialEquations: An International Journal , 23(6):1343–1367, 2007.3122] T. Xiong, J.-M. Qiu, and Z. Xu. A parametrized maximum principle preserving flux limiterfor finite difference RK-WENO schemes with applications in incompressible flows.
Journal ofComputational Physics , 252:310–331, 2013.[23] T. Xiong, J.-M. Qiu, and Z. Xu. High order maximum principle preserving discontinuousGalerkin method for convection-diffusion equations.
SIAM Journal on Scientific Computing ,37(2):583–608, 2015.[24] T. Xiong, J.-M. Qiu, Z. Xu, and A. Christlieb. High order maximum principle preserving semi-Lagrangian finite difference WENO schemes for the Vlasov equation.
Journal of ComputationalPhysics , 273:618–639, 2014.[25] H. Yoshida. Construction of higher order symplectic integrators.
Physics letters A , 150(5-7):262–268, 1990.[26] X. Zhang and C.-W. Shu. On maximum-principle-satisfying high order schemes for scalarconservation laws.
Journal of Computational Physics , 229(9):3091–3120, 2010.[27] X. Zhang and C.-W. Shu. On positivity preserving high order discontinuous Galerkin schemesfor compressible Euler equations on rectangular meshes.
Journal of Computational Physics ,229:8918–8934, 2010.[28] L. Zhou, Y. Xia, and C.-W. Shu. Stability analysis and error estimates of arbitrary LagrangianEulerian discontinuous Galerkin method coupled with Runge Kutta time-marching for linearconservation laws.