{'fraction_non_alphanumeric': 0.05418928487137125, 'fraction_numerical': 0.023176776020769414, 'mean_word_length': 4.0757067561092475, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A mass conservative Eulerian-Lagrangian Runge-Kutta discontinuous Galerkin method for wave equations with large time stepping

We propose an Eulerian-Lagrangian (EL) Runge-Kutta (RK) discontinuous Galerkin (DG) method for wave equations. The method is designed based on the ELDG method for transport problems [J. Comput. Phy. 446: 110632, 2021.], which tracks solution along approximations to characteristics in the DG framework, allowing extra large time stepping sizes with stability. The wave equation can be written as a first order hyperbolic system. Considering each characteristic family, a straightforward application of ELDG will be to transform to the characteristic variables, evolve them on associated characteristic related space-time regions, and transform them back to the original variables. However, the mass conservation could not be guaranteed in a general setting. In this paper, we formulate a mass conservative semi-discrete ELDG method by decomposing each variable into two parts, each of them associated with a different characteristic family. As a result, four different quantities are evolved in EL fashion and recombined to update the solution. The fully discrete scheme is formulated by using method-of-lines RK methods, with intermediate RK solutions updated on the background mesh. Numerical results on 1D and 2D wave equations are presented to demonstrate the performance of the proposed ELDG method. These include the high order spatial and temporal accuracy, stability with extra large time stepping size, and mass conservative property.

One effective numerical approach for solving the hyperbolic equation is the characteristic method, which updates time-dependent solution by tracking characteristics. For scalar hyperbolic equations, there have been many pioneering works. In [7], Celia etc. developed an Eulerian Lagrangian Localized Adjoint Methods (ELLAM) [7], which introduces an adjoint problem for each test function in the continuous finite element framework and has been applied to different problems [34,33]. On the other hand, the ELDG [35,6] and SLDG [4] are being developed in the discontinuous Galerkin finite element framework with a similar introduction of adjoint problems for test functions. Another line of development, that is closely related to this work, is the Arbitrary Lagrangian Eulerian (ALE) DG method [20,17]. Both ELDG and ALE-DG 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 is designed for better shock resolution.
In this paper, we propose a mass conservative EL RK DG method for hyperbolic systems (1.2) with large time stepping sizes. We start from 1D cases, for which we consider characteristic variables and the associated characteristic space-time regions. For each characteristic family, a straightforward application of ELDG will be to transform to characteristic variables, evolve them on associated space-time regions, and transform them back to original variables. However, the mass conservation could not be guaranteed in a general setting. We decompose each variable into two parts, each of them associated with different characteristic families; as a result, four different quantities are evolved in EL fashion and recombined to update the solution. The fully discrete scheme is formulated by using method-of-lines RK methods, with intermediate RK solutions updated on the background mesh. For 2D hyperbolic systems, characteristic variables are no longer constant along the bicharacteristics, characteristic Galerkin [30] or evolution Galerkin [3,22,23] methods have been proposed by taking into account information propagated in all bicharacteristic directions. However, the algorithm implementation is very complex. In this paper, we use the splitting method for higher dimensional problem, which maintain the simplicity of ELDG method for 1D cases, as well as other great properties such as mass conservation, high order spatial and temporal accuracy, and allows for extra large time steps with stability.
This paper is organized as follows. In Section 2, we review the ELDG for one-dimensional (1D) linear transport problems. In Section 3, we develop the ELDG method for a first-order hyperbolic system by evolving each component associated with each characteristic families, recombining them to update the solution. We also develop it to 2D problems by dimensional splitting. Mass conservation of ELDG schemes are proved. In Section 4, performance of the proposed ELDG method is shown through extensive numerical tests. Finally, concluding remarks are made in Section 5.

Review of ELDG formulation for 1D linear transport problems
To illustrate the key idea of the ELDG scheme, we start from a 1D linear transport equation in the following form For simplicity, we assume periodic boundary conditions, and the velocity field a(x, t) is a continuous function of space and time. We perform a partition of the computational domain x a = x 1 and define ∆x = max j ∆x j . We define the finite dimensional denotes the set of polynomials of degree at most k. We let t n be the n-th time level and ∆t = t n+1 − t n to be the time-stepping size.
The key idea in the ELDG formulation is, design adjoint problems for test functions to take advantage of information propagation along characteristics. The ELDG method proposed in [6] is being the dynamic interval, see , t n+1 ) and I j .
] is the upstream cell of I j at t n . The dynamic interval ofĨ j (t) can always be linearly mapped to a reference cell ξ in I j by a mapping. A local adjoint problem of (2.1) for the test function is defined as: Notice that the test function ψ stays the same polynomial, ifĨ j (t) is mapped to a reference interval where F (u) . = (a − α)u andF is the Lax-Friedrichs flux. A method-of-lines RK discretization can be used for high order temporal accuracy [6].
3 The ELDG algorithm for wave equation.
The design of the ELDG algorithm for wave equation shares a similar spirit as the 1D scalar case.
We start from a 1D wave equation. We first rewrite it as a linear system and formulate the ELDG scheme by tracking information along different characteristics families.

1D wave equation
We consider the wave equation For simplicity, we assume periodic boundary conditions, and the velocity field a(x) is a continuous and periodic function of space. Let u 1 = u t and u 2 = u x , then we can rewrite (3.1) as a linear We use the following notations for the eigen-decomposition of A(x): contains the right eigenvectors, and contains the left eigenvectors.
In the following, we propose a mass conservative ELDG scheme for the system (3.2) by the procedure below. The notion of its background meshes is the same as the 1D scalar case.
(1) Two partitions of space-time regions Ω (1) j and Ω (2) j . According to the first and second characteristic families, we partition the computational domain as two sets of space-time regions          (2) Adjoint Problems. We consider an adjoint problem for the first characteristic family on Ω (1) j : Similarly on Ω where (3) Formulation of a mass-conservative semi-discrete ELDG scheme. Before we present a mass-conservative semi-discrete ELDG scheme, we formulate the scheme by a localized characteristic field. In particular, a piecewise constant a j approximating a(x) in (3.4) is defined on I j , and the corresponding and (3.10) Define l (1) j is locally defined on Ω (1) j approximating l (1) (x). For simplicity, we only present the first order ELDG scheme. Take the vector product of l (1) j from left with (3.2), we have a scalar equation (3.11) Next, integrating over the space-time intervalĨ (1) j (t), then we have (3.12) The first order ELDG discretization of eq. (3.12) is to find l Here " f 11 at a cell boundary can be taken as a monotone flux, e.g. the Lax-Friedrichs flux |}. Similarly, we can easily update l (2) j U h related to λ (2) in the following: where l (2) j is a constant vector, locally defined onĨ (2) j (t) approximating l (2) (x). A simple first order ELDG scheme is composed by two evolution equations (3.13) and (3.14). That is, we can update u 1 h by (3.13),(3.14) and j (t))dt, (3.15) where U n h and U n+1 h are defined on the background mesh I j , U h (t) and U (2) h (t) are defined on the space-time dynamic meshesĨ (1) j (t) andĨ (2) j (t) respectively. Similarly, we can update u 2 h . We apply forward-Euler method for time discretization with above ELDG scheme (3.15): Remark 3.1. The above ELDG scheme is not mass conservative for two reasons: (1) Flux terms can't cancel each other as r 11 j l (1) j and r 12 j l (2) j are discontinuous across cell boundary of Ω To obtain the mass conservation, a critical point is that eigen-decomposition has to be consistent among two characteristic families and independent of partitions Ω j . We choose a p (x) as a continuous approximation of a(x) which leads that R p (x), Λ p (x), R −1 p (x) are the approximations of and The following equalities hold They are critical to design a mass-conservative ELDG scheme. Take the vector product of r p 11 (x)l (1) p (x) from left with (3.2), we have a scalar equation Multiply ψ (1) to the equation above, Next, sum equations (3.21) and (3.22), and integrate over the space-time interval Ω  .23) gives Similarly, we have an equation related to λ (2) d dt Ĩ (1)  .7): As the ELDG method for scalar, the dynamic interval of I  j (t) can always be linearly mapped to a reference cell ξ ∈ I j , the ELDG discretization of eq. (3.27) is to find j (t))dt, (3.28) for ψ (1) (x, t) satisfying the adjoint problem (3.5) with ∀Ψ(x) = ψ(x, t n+1 ) ∈ P k (I j ). Here (t), t) and " f 11 at a cell boundary can be taken as a monotone flux, e.g.
the Lax-Friedrichs flux ' f 11 (t), t) and " f 12 can be similarly defined at a cell boundary. We can similarly obtain the ELDG scheme for u 2 h .
To update (3.28) from U n h to U n+1 h , we first apply the forward Euler time discretization to get 1st order accuracy, then we generalize the scheme to general RK methods. There are two main steps involved here. In order to describe the implementation procedure of the fully discrete ELDG scheme, we define the L 2 projection.
The evaluation of the right hand side of (3.30) can be done in a subinterval-by-subinterval fashion. The implementation details can be found in [16].
Then, we propose a fully discrete EL RK DG scheme with procedure as described: 1. Obtain the initial condition on upstream meshesĨ (1) j (t n ) andĨ (2) j (t n ) of (3.28) by U     (b1) For forward-Euler method, compute We compute the four integration terms of p ) x AU n h ψ (2) (x, t n )dx highlighted in blue with U n h on background meshes in order to makes the scheme satisfy mass conservation. Thus they have to be evaluated subinterval-by-subinterval since DG solution is discontinuous. In summary, we can get where t (1) = t n+1 , u 1 h (t (1) ) and u 2 h (t (1) ) are defined on background mesh I j .
(b3) For general RK methods with intermediate stages, we will update intermediate RK solutions on background mesh as in [11]. For example, for a 2nd order RK method with Butcher tableau 3.2. It has an intermediate stage at t (1) = t n+ ∆t 2 , we propose the following steps as [11], see Figure 3.3.
i. We denote the dynamic domain tracking I j from t (1) to t n with speed ν j,(1) (t) related to second characteristic is defined similarly. Then we can update U (1) h onĨ (1) j,(1) (t) andĨ (2) j,(1) (t) from t n to t (1) as in a forward-Euler method.
ii. We update U n+1 h on dynamic domainĨ  and U h (t (1) ) with projection ontoĨ (1) j,(1) (t (1) ). conservative. In particular, given a DG solution u h (x, t n ) ∈ V k h with a periodic boundary condition, we have Proof. We firstly consider the forward Euler time discretization. Taking Ψ = 1 and F = 0 in the scheme (3.31), we have

2D linear wave equation
The solution for high-dimensional hyperbolic systems is given by means of a characteristic cone, rather than individual characteristic lines. Numerically, characteristic Galerkin [30]  introduced, which is in a predictor-corrector plus finite volume framework to get higher accuracy.
FVEG method has been developed in [28] and widely applied in [21,28,24], and the stability and accuracy have been investigated in [25,28,27,26]. Even though the FVEG method can achieve high-order accuracy and stability with extra large step, the algorithm implementation is still very complex for high-dimensional problems. In this paper, we use the splitting method for higher dimensional problem.
Consider a linear 2D wave equation It can be written as a first order linear hyperbolic system

The domain Ω is partitioned into rectangular meshes with each computational cell Ω
, where we use the piecewise Q k tensor-product polynomial spaces. Then we extend ELDG algorithm to 2D problems via dimensional splitting [31]. 1. We first locate (k + 1) 2 tensor-product Gaussian nodes on cell Ω ij : (x i,p , y j,q ), p, q = 0, ..., k.
For example, see • Evolve 1D equation (3.36) at different y j,q s for a half time-step ∆t/2, see Figure 3.4 (middle). For each y j,q , the (k + 1) point values are mapped to a P k polynomial per cell, then the 1D system (3.36) is evolved by the proposed ELDG scheme. Finally, we can map the evolved P k polynomial back to the (k + 1) point values to update the solution.
• Evolve 1D system (3.37) at different x i,p s for a full time-step ∆t as above, see Figure   3.4 (right).
• Evolve 1D system (3.36) at different y j,q s for another half time-step ∆t/2.
The splitting 2D ELDG formulation maintains high order accuracy in space, extra large time stepping size with stability and mass conservation; and has a second order splitting error.

Numerical results
In this section, we show numerical results of the proposed scheme for wave equations. We set the time stepping size as ∆t = CF L a ∆x for 1D and ∆t = CF L a ∆x + b ∆y , where a and b are maximum eigenvalues of coefficient matrixes in x-and y-directions respectively. We use RK4 [18] for time discretization. We study the following aspects: the spatial order of convergence by using small enough time stepping sizes, the spatial super-convergence of the post-processed solutions [2,9] produced by convolving the ELDG solution with a suitable kernel consisting of B-splines for the purpose of improving spatial convergence rate, the temporal order of convergence and numerical stability under a large time stepping size by varying CFL for a fixed spatial mesh.    [2,9]. We can observe the optimal convergence rate k + 1 and 2k + 1. We vary time stepping sizes, with fixed well-resolved spatial meshes, and plot error vs. CF L in Figure 4.5 and 4.6 for ELDG, ELDG1, ELDG2 and ELDG3 schemes without and with post-processed technique respectively, after a long time T = 100. The plots from post-processed ELDG schemes better show the fourth order temporal convergence. ELDG2 and ELDG3 perform comparably; they have a more restricted time step constraint than ELDG1. It indicates that, stability is affected by approximations of characteristic via the space-time partition and approximation of eigenvectors. We also note that, in both Figure 4.5 and 4.6, the CFL allowed with stability is much larger than that of the RK DG method which is 1 2k+1 . Further, we verify the mass conservative property of the ELDG schemes are around machine precision and the non-mass conservative property of the NMC ELDG schemes is presented in Figure 4.7.  Table 4.4: 1D wave equation with constant coefficient. u tt = u xx with initial condition u(x, 0) = sin(x). T = 1. We use CF L = 0.3 and CF L = 0.18 with RK4 time discretization for all P 1 and P 2 respectively. The error for only u 1 = u t was shown in this table.

1D wave equations
Mesh L 1 error Order L 1 error Order L 1 error Order L 1 error Order  Table 4.5: 1D wave equation with constant coefficient. u tt = u xx with initial condition u(x, 0) = sin(x). T = 1. We use CF L = 0.3 and CF L = 0.18 with RK4 time discretization and postprocessed technique for all P 1 and P 2 respectively. The error with post-processed technique for only u 1 = u t was shown in this table.
Mesh L 1 error Order L 1 error Order L 1 error Order L 1 error Order    improved by the mesh refinement and higher order spatial approximation.  The exact solutions u 1 and u 2 are discontinuous piecewise constants with moving discontinuities.
It is a challenging test for controlling oscillations around discontinuities. We adopt a simple TVD limiter on background mesh at each RK stages with M = 0 in [10] for all schemes. As shown in with initial condition u(x, 0) = cos(x). T = 1. We use CF L = 0.1 for P 1 and P 2 with RK4 time discretization. The error for only u 1 = u t was shown in this The expected optimal spatial accuracies of the ELDG methods without and with post-processing technique are shown in Table 4.6 and Table 4.7 respectively. In Figure 4.11 and 4.12, we plot the L ∞ error versus CF L of ELDG methods without and with post-processing technique respectively.
The following observations are made: (1) The high order accuracy of the RK method reduce the error magnitude when large time stepping size is used; (2) The ELDG methods with RK4 time discretization perform well around and before CF L = 1, which is well above the stability constraint of the RK DG method 1/(2k + 1) for P k approximations. (3) After CF L = 1 and before stability constraint of the method, the temporal convergence order is observed to be consistent with the Table 4.7: 1D wave equation with variable coefficient and source term. u tt = (a 2 (x)u x ) x + f (x, t) with initial condition u(x, 0) = cos(x). T = 1. We use CF L = 0.1 for P 1 and P 2 with postprocessed technique and RK4. The error for only u 1 = u t was shown in this    order wave equation u tt = u xx + u yy can be written as the following first order linear hyperbolic system: with period boundary conditions in both x and y directions. The exact solution is We notice that the two matrices in equation (4.2) don't commute, thus the linear system can not be reduced to 2-D scalar problems. We test accuracy for Q k ELDG methods with RK4 and 4th splitting method [5] at T = 1 for k = 1, 2 with CF L = 0.1 in Table 4.8. As expected, the (k + 1)th order convergence is observed for these methods. We plot the L ∞ error versus CF L of ELDG methods with Q 1 (left) and Q 2 (right) polynomial spaces for this case with Strang splitting and 4th order splitting in Figure 4.13, which shows that second and forth order splitting errors are dominant when time-stepping sizes are large enough. The CFL constraint with stability for ELDG method is larger than that for general RK DG method when high order time discretization is applied.   . T = 1.
linear hyperbolic system: The simplest case is to choose a(  For this example, we show the numerical ELDG Q 2 solution u 1 at times T = 0.5; 1; 1.5; 2 in Figure   4.14. Next we consider the system (4.4) with the initial condition      u 1 (x, y, 0) = 2 cos(x + y), u 2 (x, y, 0) = cos(x + y), u 3 (x, y, 0) = cos(x + y), where a(x, y) = 1+0.5 sin(x+y), b(x, y) = (4 − (1 + 0.5 sin(x + y)) 2 ) and the boundary condition is periodic in both x and y directions. The exact solution is      u 1 (x, y, t) = 2 cos(x + y + 2t), u 2 (x, y, t) = cos(x + y + 2t), u 3 (x, y, t) = cos(x + y + 2t). We report the spatial accuracy of Q k ELDG methods in Table 4.9. The expected optimal convergence is observed. We plot the L ∞ error versus CF L of ELDG methods in Figure 4.15. The ELDG Figure 4.14: Plots of the ELDG numerical solutions u 1 and their contour plots at T = 0.5; 1; 1.5; 2 for 2D wave equation u tt = u xx + u yy with Gaussian function initial condition. The mesh of 80 × 80 is used with 4th splitting method and RK4 time discretization, CF L = 1. methods perform as well as that for the linear system with constant coefficient matrices, and the CFL allowed with stability is much larger than that of the RK DG method. Table 4.9: Two-dimensional linear system with variable matrices. Q k ELDG methods (k = 1, 2) with RK4 and 4th splitting time discretization methods for (4.4) with the smooth initial condition at T = 0.1, CF L = 0.1.

Conclusion
In this paper, we have developed a mass-conservative Eulerian-Lagrangian discontinuous Galerkin Future works include further theoretic development and application of ELDG methods for nonlinear hyperbolic problems.

Appendix (The notation of test function)
We give some notations which is used in the implementation of the fully discrete ELDG scheme with RK time discretization. For convenience, we only give the definitions related to Ω (1) j below because we can similarly get the definitions related to Ω (2) j . As scalar case, we take test function ψ (1) (x, t) as ψ (1) j,m (x, t) = Ψ j,m (x − α (1) (t − t n+1 )) in the adjoint problem which is a set of basis of P k (Ĩ (1) j (t)), whereĨ (1) j (t) is donated by a domain related λ (1) asĨ j (t) in scalar problem. Here, we also take Ψ j,m as orthogonal basis on I j , and let j,0 (x, t)dx, ...,