Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for conservation laws on moving simplex meshes
Pei Fu Gero Schnücke Yinhua Xia





Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for  conservation laws on moving simplex meshes

    Pei Fu[ School of Mathematical Sciences, University of Science and Technology of China, Hefei, Anhui 230026, P.R. China.
  Email: sxfp2013@mail.ustc.edu.cn. ],
Gero Schnücke[ University of Cologne, Weyertal 86-90, 50931 Köln, Email: gschnuec@math.uni-koeln.de],
 Yinhua Xia[Corresponding author.
School of Mathematical Sciences, University of Science and Technology of China,
  Hefei, Anhui 230026, P.R. China.
  Email: yhxia@ustc.edu.cn. Research supported by NSFC grants No. 11871449 and No. 11471306, and a grant from the Science & Technology on Reliability & Environmental Engineering Laboratory (No. 6142A0502020817).
	] 
     
======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




In Klingenberg, Schnücke and Xia (Math. Comp. 86 (2017), 1203-1232) an arbitrary Lagrangian-Eulerian discontinuous Galerkin (ALE-DG) method to solve conservation laws has been developed and analyzed. In this paper, the ALE-DG method will be extended to several dimensions. The method will be designed for simplex meshes. This will ensure that the method satisfies the geometric conservation law, if the accuracy of the time integrator is not less than the value of the spatial dimension. For the semi-discrete method the L^2-stability will be proven. Furthermore, an error estimate which provides the suboptimal (k+1/2) convergence with respect to the L^∞(0,T;L^2(Ω))-norm will be presented, when an arbitrary monotone flux is used and for each cell the approximating functions are given by polynomials of degree k. The two dimensional fully-discrete explicit method will be combined with the bound preserving limiter developed by Zhang, Xia and Shu in (J. Sci. Comput. 50 (2012), 29-62). This limiter does not affect the high order accuracy of a numerical method. Then, for the ALE-DG method revised by the limiter the validity of a discrete maximum principle will be proven. The numerical stability, robustness and accuracy of the method will be shown by a variety of two dimensional computational experiments on moving triangular meshes.
Key Words: Arbitrary Lagrangian-Eulerian discontinuous Galerkin method, conservation laws, moving simplex meshes, geometric conservation law, L^2-stability, error estimates, maximum principle.





§ INTRODUCTION

The present paper investigates the development and analysis of an arbitrary Lagrangian-Eulerian discontinuous Galerkin (ALE-DG) method for scalar conservation laws in several space dimensions

    ∂_tu+∇·f(u)=0,  in Ω×(0,T),

with initial condition u_0( x) and suitable boundary conditions. The domain Ω⊆^d is an open convex polyhedron and the flux f(u):=(f_1(u),...,f_d(u))^T is a suitable vector field. In general the problem (<ref>) has no classical solutions. Discontinuities like shock waves could appear in the solution, regardless of the smoothness of the initial data. Hence, the problem needs to be investigated with a class of generalized solutions. The existence of an unique physical relevant solution for the problem was proven by Kružkov in <cit.>. This solution is called entropy solution. In particular, Kružkov proved that the unique entropy solution satisfies the maximum principle. This means the entropy solution is bounded by the interval [m,M], where

    M:= x∈Ωmax u_0( x)  and   m:= x∈Ωmin u_0( x).


The discontinuous Galerkin (DG) method, introduced by Reed and Hill <cit.> in the context of a neutron transport equation, is a finite element method with discontinuous basis functions. The choice of discontinuous basis functions gives the method a local structure (elements only communicate with immediate neighbors) and the property to handle complex mesh geometries. These features of the DG method are attractive for parallel and high performance computing. Therefore, in particular, the explicit Runge-Kutta DG (RK-DG) method for convection-dominated problems developed and analyzed by Cockburn, Shu and several co-authors in a series of publications (cf. the review article <cit.> for a summation of their pioneering works) became very popular in the last decades.

The RK-DG method of Cockburn and Shu was developed for a static computational mesh, but in engineering applications like aeroelastic computations of wings (cf. for instance Robinson et al. <cit.>) numerical methods with a deformable moving mesh are desirable. Nevertheless, a deformable computational domain can lead to strong distortions in the mesh geometry which can be the source of numerical artifacts and instabilities. In the arbitrary Lagrangian-Eulerian (ALE) approach the mesh can move with the fluid like in the Lagrangian specification or the mesh can be static as in the Eulerian specification. This flexibility has a stabilizing effect on an ALE method, since it is possible to switch to the Eulerian specification whenever distortions appear in the mesh geometry. The ALE kinematics were rigorously described by Donea et al. <cit.>. Moreover, in the literature there are different strategies to combine the ALE approach with the RK-DG method. Among others Lomtev, Kirby, Karniadakis <cit.>, Nguyen <cit.>, Persson et al. <cit.>, Kopriva et al. <cit.> and Boscheri, Dumbser <cit.> developed and analyzed ALE-DG methods for convection-dominated problems on a moving domain.

In this paper, an ALE-DG method for solving the problem (<ref>) on moving simplex meshes is introduced. This method is an extension of the ALE-DG method developed by Klingenberg et al. <cit.>. In order to describe the ALE kinematics, we assume that the distribution of the grid points is explicitly given for an upcoming time level by a suitable moving grid methodology. On the basis of this assumption, we can define local affine linear ALE mappings which connect the time-dependent simplex cells with a time-independent reference simplex cell. This simple construction of the ALE mappings ensures that our ALE-DG method has a local structure like the RK-DG method and the discrete geometric conservation law (D-GCL) is satisfied, when a suitable high order accurate Runge-Kutta (RK) method is used. The geometric conservation law (GCL) describes the time evolution of the metric terms in a grid deformation method and has an important influence on the stability and accuracy of a method. The significance of the GCL was first analyzed by Lombard and Thomas <cit.>, thenceforth the GCL was investigated in the context of moving mesh finite volume and finite element methods by Farhat et al. <cit.>, Mavriplis, Yang <cit.> and Étienne, Garon, Pelletier <cit.>.

Besides the D-GCL, a discrete maximum principle is discussed for our ALE-DG method. In general, even on a static mesh, it is not easy to design a high order method which satisfies a discrete maximum principle for the problem (<ref>) without affecting the high order accuracy of the method. Two approaches are commonly used in the literature. The first approach is the flux correction approach. Based on this approach Xu <cit.> developed a technique to ensure that a high order method satisfies the maximum principle and maintains the high order accuracy. Moreover, an algebraic flux correction approach for finite element methods was introduced by Kuzmin in <cit.> and <cit.>. Another approach was developed by Zhang and Shu <cit.>. This approach based on a bound preserving limiter which does not affect the high order accuracy of a high order method. In particular, the bound preserving limiter was developed for rectangular meshes by Zhang, Shu <cit.> and for triangular meshes by Zhang, Xia, Shu <cit.>. This approach allows the development of high order accurate maximum principle satisfying schemes by a simple investigation of the forward Euler step, since the common convexity argument (cf. Gottlieb and Shu <cit.>) can be used to extend the result for the forward Euler step to the high order total-variation-diminishing RK (TVD-RK) methods. However, Farhat, Geuzaine and Grandmont <cit.> proved that for ALE finite volume methods a discrete maximum principle is satisfied, if and only if the D-GCL is satisfied. Unfortunately, for our ALE-DG method the D-GCL is only fulfillment, if the accuracy of the RK method corresponds with the spatial dimension. Hence, we cannot expect that our forward Euler ALE-DG method satisfies a discrete maximum principle, when the bound preserving limiter is applied. Nevertheless, it turns out that the GCL is an ODE in our ALE-DG method. This ODE can be solved exactly by a RK method with an order not less than the value of the spatial dimension. In two dimensions, we use the second and the third order TVD-RK methods developed by Shu in <cit.> to solve the GCL and the actual ALE-DG method. Then the RK stage solutions for the GCL  are used to update the metric terms in the RK stages of the actual ALE-DG method. This time integration strategy allows to develop second and third order accurate fully-discrete ALE-DG methods. We prove that these methods satisfy a discrete maximum principle when the bound preserving limiter is applied. Furthermore, we present numerical experiments which support the expectation that the ALE-DG method also satisfies a discrete maximum principle when the five stage fourth order TVD-RK method developed by Spiteri and Ruuth in <cit.> and the bound preserving limiter are used.

In addition, we present an a priori error estimate for our ALE-DG method. A priori error analysis to smooth solutions of the second and third order RK-DG method on a static mesh for scalar and symmetrizable systems of conservation laws were mainly done by Zhang, Shu et al. More precisely, Zhang and Shu proved in <cit.> and <cit.> that under a slightly more restrictive Courant-Friedrichs-Lewy (CFL) constraint than the commonly used constraint the a priori error of the second order RK-DG method behaves as 𝒪( t^2+h^k+1/2) in the L^2-norm, when a local polynomial basis of degree k≥ 1 and an arbitrary monotone flux are applied. In this context the quantity t denotes the time step and h denotes the maximum cell length. Likewise, Zhang and Shu proved in <cit.> that under the usual CFL constraint the a priori error of the third order RK-DG method for scalar conservation laws behaves as 𝒪( t^3+h^k+1/2) in the L^2-norm. This result was extended to symmetrizable systems of conservation laws by Luo, Shu, and Zhang <cit.>. Furthermore, a priori error estimates for the third order RK-DG method in the context of linear scalar conservation Laws with discontinuous initial data were proven in <cit.>. Error estimates for fully discrete ALE-DG methods to solve linear conservation laws were proven by Zhou, Xia and Shu in <cit.>. In this work, we merely prove that for smooth solutions of the problem (<ref>) the a priori error of the semi-discrete ALE-DG method behaves as 𝒪(h^k+1/2), when polynomials of degree k≥max{ 1,d/2} are used on the reference cell and an arbitrary monotone flux is applied.

The rest of the paper is organized as follows. In Section 2, we introduce the local affine linear ALE mappings, a time-dependent test function space and our semi-discrete ALE-DG method. In Section 3, we present some theoretical results for the semi-discrete ALE-DG method. In particular, the L^2-stability is proven. Afterward, in Section 4, the fully-discrete ALE-DG method is investigated. We prove that the D-GCL is satisfied under certain conditions which are related to the spatial dimension. Furthermore, in two dimensions second and third order accurate fully-discrete ALE-DG methods on moving triangular meshes are presented. We prove that these methods satisfy a discrete maximum principle, when the bound preserving limiter for triangular meshes developed by Zhang, Xia, Shu in <cit.> is applied. In Section 5, we validate the theoretical results by some computational examples and show that the ALE-DG method is numerically stable and high order accurate. Finally, we give some concluding remarks in Section 6.



  §.§.§ Constants and notation

In the present paper, vectors, vector valued functions and matrices are denoted by bold letters. Scalar quantities are denoted by regular letters. The set K(t) denotes a time-dependent open simplex cell in a d dimensional domain with the edges F_K(t)^ν, ν=1,…,d+1. Volume integrals with respect to the open set K(t) and surface integrals with respect to the edges F_K(t)^ν, ν=1,…,d+1, are denoted by the bracket notation. Hence, for all v,w∈L^2(K(t))∪L^2(∂ K(t)) and ν=1,…,d+1 the notations (v,w)_K(t):=∫_K(t)vw   dx, ⟨ v,w⟩ _F_K(t)^ν:=∫_F_K(t)^νvw dΓ and ⟨ v,w⟩ _∂ K(t):=∑_ν=1^d+1⟨ v,w⟩ _F_K(t)^ν are applied. Furthermore, to avoid confusion with different constants, we denote by C a positive constant, which is independent of the mesh size and the numerical solutions for the conservation law (<ref>), but it may depend on the solution of the PDE and may have a different value in each occurrence.




§ THE ALE-DG DISCRETIZATION

In this section, we present the semi-discrete ALE-DG discretization of the problem (<ref>). At first, the ALE framework to derive the ALE-DG method in several dimensions is briefly listed. Afterward, the ALE framework is used to derive the semi-discrete ALE-DG method for solving the problem (<ref>).



 §.§ The ALE-DG setting

In this section, we present the time dependent cells and introduce some identities for the metric quantities to transform derivatives on a reference cell.



  §.§.§ The time-dependent simplex mesh

We assume that there exists a regular mesh 𝒯_(t_n) of simplices at any time level t_n, n =0,…,𝒩, which covers exactly the convex polyhedron domain Ω such that

    Ω=⋃{K(t_n)| K(t_n)∈𝒯_(t_n)}.

The mesh topology of 𝒯_(t_n) and 𝒯_(t_n+1) is assumed to be the same. This means:


  a) 𝒯_(t_n) and 𝒯_(t_n+1) are simplex meshes of the domain Ω.

  b) 𝒯_(t_n) and 𝒯_(t_n+1) have the same number of cells.

  c) The cells of both simplex meshes are positively oriented with respect to the reference simplex

    :={ξ=(ξ_1,...,ξ_d)^T∈^d: ξ_ν≥0, ∀ν, and ∑_ν=1^dξ_ν≤1}.

		
The d+1 vertices of each simplex K(t_n)∈𝒯_(t_n) are denoted by v_1^n,…, v_d+1^n. We define for all t∈[t_n,t_n+1] and ℓ =1,...,d+1 time-dependent straight lines

    v_ℓ(t):= v_ℓ^n+ω_K^n,ℓ(t-t_n),    ω_K^n,ℓ:=1/ t( v_ℓ^n+1- v_ℓ^n).

These straight lines are for any t∈[t_n,t_n+1] the vertices of a time-dependent simplex cell given by

    K(t):=   int(conv{ v_1(t),…, v_d+1(t)}),
      ∂ K(t)=⋃_ν=1^d+1F_K(t)^ν,  
        F_K(t)^ν:=conv({ v_1(t),…, v_d+1(t)}∖{ v_ν(t)}),

where int(·) and conv(·) denote the interior and the convex hull of a set. In the following, the set of all time-dependent cells K(t) is denoted by 𝒯_(t).
Furthermore, for any cell K(t)∈𝒯_(t) the diameter of the cell and the radius of the largest ball, contained in K(t), are denoted by h_K(t) as well as ρ_K(t). Additionally, we define the following global length

    h:=t∈[0,T]maxK(t)∈𝒯_(t)maxh_K(t).

Henceforth, we assume:


  (A1) The domain Ω is for all  t∈[0,T] exactly covered by the time-dependent cells (<ref>) such that Ω=⋃_K(t)∈𝒯_(t)K(t).

  (A2) For all  t∈[0,T] and all cells K(t)∈𝒯_(t) J_K(t)=det( A_K(t))>0.

  (A3) It exists constants κ>0 and τ>0, independent of h, such that for all  t∈[0,T]
 
    h_K(t)≤κρ_K(t)  and   h ≤τ h_K(t),  ∀ K(t)∈𝒯_(t).





  §.§.§ The ALE mapping and the grid velocity field

The time-dependent simplex cells (<ref>) can be mapped to the time-independent reference simplex element (<ref>) by the affine linear time-dependent mapping

    χ_K(t):→K(t),  ξ↦χ_K(t)(ξ,t):= A_K(t)ξ+ v_1(t),

where the matrix A_K(t) is given by

    A_K(t):=( v_2(t)- v_1(t),…, v_d+1(t)- v_1(t)).

We note that the matrix A_K(t) is the Jacobian matrix of the mapping χ_K(t) and the corresponding determinant is

    J_K(t)=det( A_K(t))=d!|K(t)|,

where |K(t)| denotes the volume of the cell K(t). In particular, J_K(t) is independent of the spatial variables and belongs to P^d([t_n-1,t_n]). It is worth to mention that in general for non-simplicial moving meshes J_K(t) depends on spatial and temporal variables, since the shape of the elements can change when the corners move with different speed. In Figure <ref> the two dimensional situation for a triangular element and a rectangular element is illustrated. The implementation of metric quantities  which depend on spatial variables is not easy and requires caution (e.g. cf. Kopriva  <cit.>).


Since the matrix A_K(t) is the Jacobian matrix of the mapping (<ref>), 
we have the following metric transformations

    ∇·f=∇_ξ·[( A_K(t)^-1)f^*],      ∇ u= A_K(t)^-T∇_ξu^*,
        _K(t)
    =J_K(t) A_K(t)^-T_,

where f:→^d is an arbitrary vector field  with f^*=f∘χ_K(t), u is a scalar function with u^*=u∘χ_K(t), _K(t) is the normal of the cell K(t) and _ is the reference normal. A proof of these metric transformations is given in Ciarlet <cit.>. Moreover, the mapping (<ref>) provides the grid velocity field in the point x=χ_K(t)(ξ,t)

    ω_K(t)( x,t):=d/dt(χ_K(t)(ξ,t)).

The definition (<ref>) provides the following relation

    ∂_ξ_iω_j(χ_K(t)(ξ,t)),t)=[d/dt( A_K(t))]_ji,     i,j=1,…,d,

where ω_j are the coefficients of the grid velocity and [d/dt( A_K(t))]_ji are the coefficients of the matrix d/dt( A_K(t)). We are also interested to find an identity for the time derivative of the determinant J_K(t). Hence, we apply Jacobi's formula (cf. Bellman <cit.>), using the equation relating the adjugate of A_K(t) to the inverse A_K(t)^-1 and apply the identity (<ref>). This results in the identity

    d/dt(J_K(t))
    =   tr[adj( A_K(t))d/dt( A_K(t))] 
    
    =   ∑_i=1^d∑_j=1^d(-1)^i+jM_K(t)^ij[d/dt( A_K(t))]_ji
    
    =   (∑_i=1^d∑_j=1^d[ A_K(t)^-1]_ij(∂_ξ_iω_j))J_K(t)
    
    =   (∑_i=1^d[ A_K(t)^-1(∂_ξ_iω)]_i)J_K(t),

where tr(·) denotes the trace of a matrix, adj(·) denotes the adjoint of a matrix and M_K(t)^ij  is the (i,j) minor of A_K(t). Moreover, since the matrix A_K(t) does not depend on spatial variables, we obtain

    d/dt(J_K(t))=(∑_i=1^d[ A_K(t)^-1(∂_ξ_iω)]_i)J_K(t)=(∇_ξ·[ A_K(t)^-1ω])J_K(t).

Finally, we summarize some properties of the grid velocity. These properties will be used in the next sections.
 The grid velocity ω defined by (<ref>) has the properties:


  (i) For all t∈[t_n,t_n+1] the grid velocity belongs to the space  P^1(,^d).

  (ii) The grid velocity is time-independent for all points contained in the set ∂.

  (iii) The divergence of the grid velocity satisfies

    (∇·ω)J_K(t)=(∇_ξ·[ A_K(t)^-1ω])J_K(t)∈ P^d-1([t_n-1,t_n]).




Since the matrix d/dt( A_K(t)) is time independent, the property (i) follows  from the definition of the grid velocity (<ref>).

The property (ii) follows directly from the definitions of the time-dependent straight lines (<ref>), the cells (<ref>) and the grid velocity (<ref>).

Finally, the property (iii) follows from the metric transformations (<ref>) and the identity (<ref>), since [d/dt( A_K(t))]_ij is time independent and the minors M_K(t)^ij, i,j=1,…,d, belong to the space P^d-1([t_n-1,t_n]).





 §.§ The approximation space

We define the approximation space

    𝒱_h(t):={ v∈L^2(Ω): v∘χ_K(t)∈ P^k(), ∀ K(t)∈𝒯_(t)} ,

where P^k() denotes the space of polynomials in  of degree at most k. The functions from the space 𝒱_h(t) are discontinuous along the interface of two adjacent cells. Thus, we define for a function v∈𝒱_h(t), an arbitrary cell K(t)∈𝒯_(t) and all ν =1,...,d+1 the following limits

    v^int_K(t)( x):=ε→0^+lim v( x-ε𝐧_K(t)^ν),
      
    v^ext_K(t)( x):=ε→0^+lim v( x+ε𝐧_K(t)^ν),
      ∀ x∈ F_K(t)^ν,

where the vector _K(t)^ν, ν=1,...,d+1, is the outward normal of the cell K(t) with respect to the simplex face F_K(t)^ν. Then, the cell average and jump of the function v along the simplex face F_K(t)^ν are defined by

    v:=1/2(v^int_K(t)+v^ext_K(t)),
      v:=v^ext_K(t)-v^int_K(t).





 §.§ The semi-discrete ALE-DG method

For each cell K(t) ∈𝒯_(t), we approximate the solution u of the problem (<ref>) by the function

    u_h( x,t)=∑_j=1^ru_j^K(t)(t)ϕ_j^K(t)( x,t),  for all t∈[t_n,t_n+1) and  x∈ K(t),

where r:=(k+d)!/d!k! and {ϕ_1^K(t)( x,t),...,ϕ_r^K(t)( x,t)} is a basis of the space 𝒱_h(t) in the cell K(t).
The coefficients u_1^K(t)(t),...,u_r^K(t)(t) in (<ref>) are the unknowns of the method. In order to determine these coefficients, we plug the function (<ref>) in (<ref>), multiply the equation by a test function v ∈𝒱_h(t) and use the change of variables theorem for integrals. This results in the equation
 
    (J_K(t)(∂_tu_h),v^*)_+(J_K(t)(∇_x·f(u_h)),v^*)_=0,

where v^*=v∘χ_K(t). The chain rule formula and the metric transformations (<ref>) provide

    d/dtu_h^*=∂_tu_h+ω·∇ u_h=
    ∂_tu_h+ω· A_K(t)^-T∇_ξu_h^*,

where u_h^*=u_h∘χ_K(t). Next the identities (<ref>) and (<ref>) provide

    J_K(t)(∂_tu_h)=d/dt(J_K(t)u_h^*)-J_K(t)∇_ξ·[ A_K(t)^-1(ωu_h^*)].

Then the equation (<ref>) becomes

    (d/dt(J_K(t)u_h^*),v^*)_+(J_K(t)(∇_ξ·g̃(ω,u_h^*)),v^*)_=0,

where g̃(ω,u_h^*):= A_K(t)^-1(g(ω,u_h^*)) with

    g(ω,u):=f(u)-ω( x,t)u.

At this point, we proceed similar as in the derivation of the standard DG method on a static mesh. First, we apply the integration by parts formula in the second integral in (<ref>). Then, we replace the flux function g̃(ω,u_h^*,int_)·_ in the surface integrals by a numerical flux function g(ω,u_h^*,int_,u_h^*,ext_,J_K(t)(t)) with (t)=A_K(t)^-T_. The numerical flux function needs to satisfy certain properties. These properties are discussed in the Section <ref>. Finally, on the reference cell, the semi-discrete ALE-DG method appears as the following problem:

[The semi-discrete ALE-DG method on the reference cell]
Find a function u_h∈𝒱_h(t), such that for all v∈𝒱_h(t) and all cells K(t)∈𝒯_(t) holds

    (d/dt(J_K(t)u_h^*),v^*)_
     =   (J_K(t)g̃(ω,u_h^*),∇_ξv^*)_
        -⟨g(ω,u_h^*,int_,u_h^*,ext_,J_K(t)(t)),v^*,int_⟩ _∂,

where (t)=A_K(t)^-T_.

Since the test functions v^*=v∘χ_K(t) are time independent on the reference cell , we obtain

    (d/dt(J_K(t)u_h^*),v^*)_=d/dt(J_K(t)u_h^*,v^*)_=d/dt(u_h,v)_K(t).

Therefore, on the physical domain,  <Ref> is equivalent to:

[The semi-discrete ALE-DG method]
Find a function u_h∈𝒱_h(t), such that for all v∈𝒱_h(t) and all cells K(t)∈𝒯_(t) holds

    d/dt(u_h,v)_K(t)=   (g(ω,u_h),∇ v)_K(t)
       -⟨g(ω,u_h^int_K(t),u_h^ext_K(t),_K(t)),v^int_K(t)⟩ _∂ K(t).






 §.§ The numerical flux function

The numerical flux in the ALE-DG method should satisfy:


  (P1) The flux is consistent with g(ω,u)·_K(t).

  (P2) The function u↦g(ω,u,·,_K(t)) is increasing and Lipschitz continuous.

  (P3) The function u↦g(ω,·,u,_K(t)) is decreasing and Lipschitz continuous.

  (P4) The flux is conservative such that

    g(ω,u_h^int_K(t),u_h^ext_K(t),_K(t))=-g(ω,u_h^ext_K(t),u_h^int_K(t),-_K(t)).



The properties (P1) - (P4) of the numerical flux supply for any cell K(t)∈𝒯_(t) and all v∈[min(u_1,u_2),max(u_1,u_2)]

    (g(ω,v)·_K(t)-g(ω,u_1,u_2,_K(t)))(u_2-u_1)≥0.

The inequality (<ref>) is the e-flux condition, which was introduced by Osher <cit.>.

In general every numerical flux with these properties can be used in  the ALE-DG method. A common example is the Lax-Friedrichs flux. This flux is given by


    g(ω,u_1,u_2,_K(t)):=g_+(ω,u_1,_K(t))-g_-(ω,u_2,_K(t)),


    g_±(ω,u,_K(t)):=1/2[λ^nu±g(ω,u)·_K(t)],


    λ^n:=max{|(∂_ug(ω,u))·_K(t)|: u∈[m,M], t∈[t_n,t_n+1]}.


We note that the functions g_±(ω,u,_K(t)) are increasing in the second argument.





§ THEORETICAL RESULTS FOR THE SEMI-DISCRETE METHOD

In this section, we present some theoretical results for the semi-discrete ALE-DG method. We start with a proof for the L^2-stability of the semi-discrete ALE-DG method. This proof requires techniques which were introduced by Jiang and Shu in <cit.> to proof a cell entropy inequality for the DG method. We note that for scalar conservation laws the function η(u)=1/2u^2 is an entropy. Thus, the L^2-stability provides also entropy stability in the sense that the total entropy at a certain time point is bounded by the total entropy at initial time.



The solution u_h of the semi-discrete ALE-DG method satisfies for any t∈ [0,T]

    u_h(t)_L^2(Ω)≤u_h(0)_L^2(Ω)

when the problem (<ref>) is considered with periodic boundary conditions.


Let K(t)∈𝒯_(t) be an arbitrary cell. We use the ALE-DG solution u_h^*=u_h∘χ_K(t) as test function in the equation (<ref>) and obtain

    (d/dt(J_K(t)u_h^*),u_h^*)_=   (J_K(t)g̃(ω,u_h^*),∇_ξu_h^*)_
       -⟨g(ω,u_h^*,int_,u_h^*,ext_,J_K(t)(t)),u_h^*,int_⟩ _∂.

We obtain by the identity (<ref>) and the change of variables theorem for integrals

    (d/dt(J_K(t)u_h^*),u_h^*)_
    =   1/2d/dt((u_h^*)^2,J_K(t))_+1/2(J_K(t)∇_ξ·[( A_K(t)^-1)ω],(u_h^*)^2)_
    
    =   1/2d/dt(u_h,u_h)_K(t)+1/2(∇·ω,u_h^2)_K(t).

Next, we define the vector valued functions

    F(u):=(∫^uf_1(u) du,...,∫^uf_d(u) du)^T  and  G(ω,u):=F(u)-1/2ωu^2.

Then, we obtain by the metric transformations (<ref>), the functions (<ref>) and the change of variables theorem for integrals














    (J_K(t)g̃(ω,u_h^*),∇_ξu_h^*)_
    -1/2(J_K(t)∇_ξ·[( A_K(t)^-1)ω],(u_h^*)^2)_
    
    =   (J_K(t)( A_K(t)^-1)f(u_h^*),∇_ξu_h^*)_-1/2(J_K(t)( A_K(t)^-1)ω,∇_ξ(u_h^*)^2)_
       -1/2(J_K(t)∇_ξ·[( A_K(t)^-1)ω],(u_h^*)^2)_
    
    =   (J_K(t)( A_K(t)^-1)f(u_h^*),∇_ξu_h^*)_-1/2(J_K(t)∇_ξ·[( A_K(t)^-1)ω(u_h^*)^2],1)_
    
    =   (J_K(t)∇_ξ·(( A_K(t)^-1)G(ω,u_h^*)),1)_=(∇·G(ω,u_h),1)_K(t)

Furthermore,  the divergence theorem gives

    (∇·G(ω,u_h),1)_K(t)
    =∑_ν=1^d+1⟨G(ω,u_h^int_K(t))·_K(t)^ν,1⟩ _F_K(t)^ν,

where the vectors _K(t)^ν, ν=1,…,d+1,  are the outward normals of the cell K(t) with respect to the simplex faces F_K(t)^ν, ν=1,…,d+1. Likewise, the metric transformations (<ref>) provide

    ⟨g(ω,u_h^*,int_,u_h^*,ext_,J_K(t)(t)),u_h^*,int_⟩ _∂
    
    =   ⟨g(ω,u_h^int_K(t),u_h^ext_K(t),_K(t)),u_h^int_K(t)⟩ _∂ K(t)
    
    =   ∑_ν=1^d+1⟨g(ω,u_h^int_K(t),u_h^ext_K(t),_K(t)^ν),u_h^int_K(t)⟩ _F_K(t)^ν.

Next, we rearrange the equation (<ref>) by applying the identities (<ref>), (<ref>), (<ref>) and (<ref>). This results in the equation

    1/2d/dt(u_h,u_h)_K(t)=   ∑_ν=1^d+1⟨G(ω,u_h^int_K(t))·_K(t)^ν,1⟩ _F_K(t)^ν
       -∑_ν=1^d+1⟨g(ω,u_h^int_K(t),u_h^ext_K(t),_K(t)^ν),u_h^int_K(t)⟩ _F_K(t)^ν.

Then, we sum the equation (<ref>) over all cells K(t)∈𝒯_(t) and obtain

    1/2d/dt‖ u_h‖ _L^2(Ω)^2 - 1/2∑_K(t)∈𝒯_(t)∑_ν=1^d+1⟨G(ω,u_h)·_K(t)^ν,1⟩ _F_K(t)^ν
    
    + 1/2∑_K(t)∈𝒯_(t)∑_ν=1^d+1⟨g(ω,u_h^int_K(t),u_h^ext_K(t),_K(t)^ν),u_h⟩ _F_K(t)^ν   =0,

since we consider the problem (<ref>) with periodic boundary conditions. The definition of the function G(ω,u) in (<ref>) yields ∂_uG(ω,u)=g(ω,u). Thus, for any cell K(t)∈𝒯_(t) the mean value theorem and the e-flux condition (<ref>) provide for all ν=1,…,d+1

    ⟨g(ω,θ_K(t)^ν)·_K(t)^ν-g(ω,u^int_K(t),u^ext_K(t),_K(t)^ν),u_h⟩ _F_K(t)^ν≥0,

where θ_K(t)^ν is a value between min(u_h^int_K(t),u_h^ext_K(t)) and max(u_h^int_K(t),u_h^ext_K(t)).
Hence, the inequality (<ref>) follows by integrating the equation (<ref>) over the interval [0,t].


Furthermore, for sufficiently smooth solutions of the initial value problem (<ref>), we have a suboptimal a priori error estimate in the sense of the L^∞(0,T;L^2(Ω))-norm for the semi-discrete ALE-DG method.

Consider the initial value problem (<ref>) with periodic boundary conditions and let u∈W^1,∞(0,T;H^k+1(Ω)) be the exact solution, the flux function f∈𝒞^3(,^d) and the grid velocity field ω be bounded in Ω×[0,T] and have bounded derivatives and u_h be the solution of the semi-discrete ALE-DG method. The test function space (<ref>) is given by piecewise polynomials of degree k≥max{ 1,d/2}. Furthermore, (A1) - (A3) are satisfied and the global length h is given by (<ref>).
Then, it exists a constant C, which depends on the final time T and is independent of h and u_h, such that 

    ‖ u-u_h‖ _L^∞(0,T;L^2(Ω))≤ Ch^k+1/2.


<Ref> can be proven with standard techniques from approximation theory (cf. Ciarlet <cit.>). In particular, a one dimensional proof is given in <cit.> and in <cit.> two dimensional error analysis for a semi-discrete ALE-DG method to solve the Hamilton-Jacobi equations is given. In addition, error analysis for the fully-discrete ALE-DG methods is given in <cit.>. Since, there are already these publications on error analysis for the ALE-DG method in the literature and the proof  of <Ref> vary merely in technical details, we skip the proof in this paper.

The proof of <Ref> requires the a priori assumption

    ‖ u-u_h‖ _L^∞(Ω×[0,T])≤ h^1/2.

These a priori assumption is not necessary, if the problem (<ref>) is considered with a linear flux function f(u)= cu, c∈^d. In the one-dimensional case, the a priori assumption

    ‖ u-u_h‖ _L^∞(0,T;L^2(Ω))≤ h

can be applied, since for d=1 the inequality (<ref>) and the inverse inequality <cit.> provide the a priori assumption (<ref>). The assumption (<ref>) was applied in <cit.>. Moreover, the statement of <ref> can be also proven, for flux functions with less smoothness. Nevertheless, this requires a more restrictive a priori assumption and more restrictive bounds for the parameter k. This assumption was applied by Klingenberg et al. in <cit.>.





§ THE FULLY-DISCRETE METHOD

In this section, we discuss the time discretization of the ALE-DG method. In the first part of this section, we prove that the fully-discrete ALE-DG method satisfies  the discrete geometric conservation law (D-GCL). Afterward, in two dimensions, we prove that the second and the third order accurate fully-discrete ALE-DG methods satisfy the maximum principle, when the bound preserving limiter developed by Zhang, Xia and Shu in <cit.> is applied.




 §.§ Total-variation-diminishing Runge-Kutta (TVD-RK) methods

We apply the high order TVD-RK methods developed by Shu in <cit.> for the time discretization, which is also known as strong stability preserving Runge-Kutta (SSP-RK) methods <cit.>. For a given ODE u_t=q(u,t), a s-stage TVD-RK method can be written as 

    u^n,0=     u^n, 
    
    u^n,i=    ∑_j=0^i-1(α_iju^n,j+ tβ_ijq(u^n,j,t_n+γ_j)),   for i=1,…,s, 
    
    u^n+1=    u^n,s,

where t_n+γ_j:=t_n+γ_j t. The coefficients of the s-stage TVD-RK method (<ref>) need to satisfy

    0≤γ_j≤1;    α_ij,β_ij≥0;    α_ij=0  ⇔  β_ij=0;    ∑_j=0^i-1α_ij=1;

for all i=1,…,s and j=0,…,s-1. In the Section <ref>, we will investigate the commonly used second and third order TVD-RK methods from Shu <cit.>. The coefficients for these methods are given in the Tables <ref> and <ref>.


tableThe coefficients for the second order TVD-RK method from Shu <cit.>.

TVD-RK2
 
γ_0=0     α_10=1     β_10=1 

γ_1=1     α_20=1/2, α_21=1/2     β_20=0, β_21=1/2
 







tableThe coefficients for the third order TVD-RK method from Shu <cit.>.

TVD-RK3
 
γ_0=0     α_10=1     β_10=1 

γ_1=1     α_20=3/4, α_21=1/4     β_20=0, β_21=1/4  

γ_2=1/2     α_31=1/3, α_31=0, α_32=2/3      β_30=0, β_31=0, β_32=2/3
 








 §.§ The discrete geometric conservation law (D-GCL)

For the time discretization of the semi-discrete ALE-DG Problem <ref>, it is convenient to introduce the notation

    𝒢(u_h^*,v^*,J_K(t),t):=   (J_K(t)g̃(ω(t),u_h^*),∇_ξv^*)_
       -⟨g(ω,u_h^*,int_,u_h^*,ext_,J_K(t)(t)),v^*,int_⟩ _∂,

where (t)=A_K(t)^-T_. At this point, we note that the grid velocity satisfies the property (ii) in <Ref>. Therefore, in order to avoid confusion with the time-dependency of the grid velocity, we will highlight the time variable if the grid velocity depends on time, otherwise the time variable will be omitted. For all v∈𝒱_h(t) with v^*=v∘χ_K(t) the equation (<ref>) can be written as

    (d/dt(J_K(t)u_h^*),v^*)_=𝒢(u_h^*,v^*,J_K(t),t).



In accordance with Guillard and Farhat <cit.>, we introduce the following definition.

A fully-discrete moving mesh method for the initial value problem (<ref>) satisfies a D-GCL, if for all c∈ and n=0,…,𝒩-1

    u_h(𝐱,t_n)=c,  for all 𝐱∈Ω    ⇒     u_h(𝐱,t_n+1)=c,  for all 𝐱∈Ω.



Let c∈ be a constant. Then for all t∈[t_n,t_n+1], K(t)∈𝒯_(t) and v∈𝒱_h(t) with v^*=v∘χ_K(t), it holds

    𝒢(c,v^*,J_K(t),t)=(J_K(t)∇_ξ·[ A_K(t)^-1(ω(t)c)],v^*)_.



Since f(c) contains merely constant coefficients and

    g̃(ω,c)= A_K(t)^-1(g(ω,c))= A_K(t)^-1(f(c))- A_K(t)^-1(ω(t)c),

the integration by parts formula provides

    (J_K(t)g̃(ω(t),c),∇_ξv^*)_
    =     (J_K(t)∇_ξ·[ A_K(t)^-1(ω(t)c)],v^*)_
        +⟨g̃(ω(t),c)·(J_K(t)_),v^*,int_⟩ _∂.

Furthermore, since (t)=A_K(t)^-T_K_ref, the property (P1) of the numerical flux provides

    g(ω,c,c,J_K(t)(t))=g̃(ω(t),c)·(J_K(t)_).

Thus, we obtain the identity (<ref>) by (<ref>) and (<ref>).

Next, we assume that u_h^*=c solves the semi-discrete ALE-DG method Problem <ref>. Then, we obtain by (<ref>) and (<ref>)

    (d/dt(J_K(t)c),v^*)_=(J_K(t)∇_ξ·[ A_K(t)^-1(ω(t)c)],v^*)_.

The equation (<ref>) and the ODE (<ref>) are equivalent, since c is an arbitrary constant, v∈𝒱_h(t) with v^*=v∘χ_K(t) is an arbitrary test function and the quantities

    J_K(t)∈ P^d([t_n,t_n+1]),     (∇_ξ·[ A_K(t)^-1(ω(t))])J_K(t)∈ P^d-1([t_n,t_n+1])

are merely time-dependent. We note that the time evolution of the  metric terms J_K(t) needs to be respected in the time discretization of the semi-discrete ALE-DG method Problem <ref>. Therefore, we discretize the ODE (<ref>) and (<ref>) simultaneously by the same TVD-RK method. The stage solutions of the TVD-RK  discretization for (<ref>) will be used to update the metric terms in the TVD-RK discretization for (<ref>).


The fully-discrete ALE-DG method:
First, the ODE (<ref>) is discretized by a s-stage TVD-RK method:


    J_K^n,0=J_K^n, 
    
        J_K^n,i=∑_j=0^i-1(α_ijJ_K^n,j+β_ij t(∇_ξ·[ A_K^n+γ_j^-1(ω^n+γ_j)])J_K^n,j),   for i=1,…,s, 
        J_K^n+1=J_K^n,s,


where K^n+γ_j:=K(t_n+γ_j t) and ω^n+γ_j:=ω(t_n+γ_j t).
The stage solutions { J_K^n,i} _i=0^s are used to update the metric terms in the TVD-RK discretization of (<ref>). The Runge-Kutta method needs to solve the ODE (<ref>) exact such that

    J_K^n+1=J_K(t_n+1)=d!|K(t_n+1)|,    ∀ K(t_n+1)∈𝒯_(t_n+1),

where 𝒯_(t_n+1) is the regular mesh of simplices which has been used in the Section <ref> to construct the time-dependent cells (<ref>). We note that in a d-dimensional space a TVD-RK method with order greater than or equal to d is necessary to compute the metric term J_K^n+1 exactly, since the right hand side of the equation (<ref>) belongs to the space P^d-1([t_n,t_n+1]).
[The fully-discrete ALE-DG method on the reference cell]
Find a function u_h∈𝒱_h(t), such that for all v∈𝒱_h(t) there holds


    (J_K^n,0u_h^n,0,*,v^*)_=     (J_K^nu_h^n,*,v^*)_, 
    (J_K^n,iu_h^n,i,*,v^*)_=     ∑_j=0^i-1α_ij(J_K^n,ju_h^n,j,*,v^*)_
        +∑_j=0^i-1β_ij t𝒢(u_h^n,j,*,v^*,J_K^n,j,t_n+γ_j),  for i=1,…,s, 
    (J_K^n+1u_h^n+1,*,v^*)_=     (J_K^n,su_h^n,s,*,v^*)_,


where the metric terms { J_K^n,i} _i=0^s are computed by (<ref>).

We note that on a static mesh the cells K(t) are time-independent. Thus, on a static mesh, the Problem <ref> corresponds to the TVD-RK DG method developed by Cockburn and Shu in <cit.>. Next, we prove that the fully-discrete ALE-DG method satisfies the following statement.

Suppose a s-stage TVD-RK method with order greater than or equal to d is used in (<ref>) and Problem <ref>, and the solutions { J_K^n,i} _i=0^s are used to compute the metric terms in Problem <ref>. Moreover, the solution at time level t_n satisfies u_h^n,*=c∈. Then it is u_h^n,i,*=c for all i=0,…,s.


Let i∈{ 0,…,s} be an arbitrary fixed index and v∈𝒱_h(t) be an arbitrary test function. We are interested to investigate the i-th Runge-Kutta stage in the Problem <ref>. Hence, we can assume that u_h^n,j,*=c for all j=0,…,s-1. Then, the equation (<ref>) in <Ref> provides

    (J_K^n,iu_h^n,i,*,v^*)_
    =     ∑_j=0^i-1α_ij(J_K^n,jc,v^*)_
       +∑_j=0^i-1β_ij t(J_K^n,j∇_ξ·[ A_K^n+γ_j^-1(ω^n+γ_jc)],v^*)_.

Next, we multiply the equation (<ref>) by cv^* and integrate the result over the reference element . This provides the identity

    (J_K^n,ic,v^*)_
    =     ∑_j=0^i-1α_ij(J_K^n,jc,v^*)_
       +∑_j=0^i-1β_ij t(J_K^n,j∇_ξ·[ A_K^n+γ_j^-1(ω^n+γ_jc)],v^*)_.

Since the metric terms { J_K^n,j} _j=0^i-1 in (<ref>) are computed by (<ref>), the equation (<ref>) can be plugged in the equation (<ref>) and it follows

    (J_K^n,iu_h^n,i,*,v^*)_=(J_K^n,ic,v^*)_.

Thus, it follows u_h^n,i,*=c, since the metric term J_K^n,i is merely time-dependent and the test function v∈𝒱_h(t) was chosen arbitrary.

<Ref> provides that constant initial data is preserved in each Runge-Kutta stage. In particular, a direct consequence of the <Ref> is the following result.

Suppose a s-stage TVD-RK method with order greater than or equal to d is used in (<ref>) and Problem <ref>. Furthermore, the solutions { J_K^n,i} _i=0^s given by (<ref>) are used to compute the metric terms in Problem <ref>. Then the fully-discrete ALE-DG method satisfies the D-GCL in the sense of the <Ref>.





 §.§ The discrete maximum principle

In this section, we investigate the fully-discrete ALE-DG method Problem <ref> in two dimensions. First, we review the bound preserving limiter which was developed by Zhang, Xia and Shu in <cit.>.


Then, we consider the second and third order accurate TVD-RK methods given in the Tables <ref> and <ref>. We prove that for these methods the fully-discrete ALE-DG method satisfies the maximum principle when the bound preserving limiter is applied.




  §.§.§ The bound preserving limiter

In the following, we briefly review the bound preserving limiter methodology for an arbitrary s-stage TVD-RK method (<ref>). Let u_h^n be the DG solution at time level t=t_n, K⊆Ω an arbitrary cell and Q_K⊆ K be a set of quadrature points. The corresponding quadrature formula needs to be exact for the integration of polynomials with degree k on the cell K and the set Q_K needs to contain all quadrature points which are necessary to evaluate the surface integrals along the edges of K in the DG method. A quadrature rule with these properties has been developed in <cit.> and will be presented in the Section <ref>. Then the cell average

    u_K^n:=1/|K|(u_h,1)_K

of the numerical solution u_h can be computed exact by the quadrature formula. The limiter methodology is applied in two steps:


  (L1) It needs to be ensured that the s-stage TVD-RK method satisfies

    u_h^n,0|_Q_K=u_h^n|_Q_K,
       for ℓ=1,…,s: 
        u_h^n,j|_Q_K∈[m,M],   ∀ j=0,...,ℓ-1,
       ⇒  u_K^n,ℓ∈[m,M],
       u_K^n+1=u_K^n,s,

where m as well as M are given by (<ref>), u_h^n,ℓ is the solution given by the ℓ-th stage of the RK-DG method and u_K^n,ℓ denotes the corresponding cell average with respect to the cell K.

  (L2) If (L1) is satisfied, the solution u_h will be revised for all ℓ=1,...,s by

    u_h^n,ℓ|_K:=Θ(u_h^n,ℓ|_K-u_K^n,ℓ)+u_K^n,ℓ,


    Θ:=min{|M-u_K^n,ℓ/M_K-u_K^n,ℓ|,|m-u_K^n,ℓ/m_K-u_K^n,ℓ|,1},

where M_K:= x∈ Q_Kmax u_h^n,ℓ( x), m_K:= x∈ Q_Kmin u_h^n,ℓ( x). Then it is ensured that u_h^n,ℓ|_K∈[m,M], for all ℓ=1,...,s.

The stability property (L1) is the maximum principle property for the cell-averages of the numerical solution u_h. We note that by an adjustment of the CFL constraint the stability property (L1) is satisfied for any high order TVD-RK method, if the forward Euler step satisfies (L1), since the TVD-RK methods are convex combinations of the forward Euler step (cf. Gottlieb and Shu <cit.>). However, the forward Euler method is merely first order accurate and thus this method is not an appropriate choice for the discretization of the ODE (<ref>) and the semi-discrete ALE-DG method Problem <ref> in two dimensions, since it does not provide that the equation (<ref>) holds. Therefore, the convexity argument in Gottlieb and Shu <cit.> cannot be used to analyze a fully-discrete ALE-DG method which has a high order TVD-RK method as time integrator and uses the bound preserving limiter. For this reason, each high order TVD-RK method needs to be investigated separately.




  §.§.§ High order time discretization methods


The step (L1) in the bound preserving limiter methodology is related to the  cell averages of the numerical solution. Thus, a scheme satisfied by the cell averages of the ALE-DG solution u_h needs to be investigated. In this section, we will present these schemes for the TVD-RK2 and the TVD-RK3 ALE-DG methods. Later, in the Section <ref> and <ref>, these schemes will be used to prove that the TVD-RK2 and the TVD-RK3 ALE-DG methods satisfy indeed the maximum principle, when the bound preserving limiter (<ref>) is applied.  


The TVD-RK2 ALE-DG method:
We note that by the identity (<ref>) the metric term J_K(t) can be also written as 2|K(t)|. Thus, for the TVD-RK2 method in Table <ref> the scheme (<ref>) becomes


    |K^n,1|=   |K^n|+ t(∇_ξ·[ A_K^n^-1(ω^n)])|K^n|, 
    |K^n+1|=   1/2|K^n|+1/2|K^n,1|+ t/2(∇_ξ·[ A_K^n+1^-1(ω^n+1)])|K^n,1|.


Next, we apply the TVD-RK2 method and the test function v^*=1 in (<ref>) and obtain


    |K^n,1|u_K^n,1^n,1
    =   |K^n|u_K^n^n + t𝒢(u_h^n,*,1,J_K^n,t_n), 
    |K^n+1|u_K^n+1^n+1
    =   1/2|K^n|u_K^n^n+1/2|K^n,1|u_K^n,1^n,1
    + t/2𝒢(u_h^n,1,*,1,J_K^n,1,t_n),


where u_K^n^n, u_K^n,1^n,1 and u_K^n+1^n+1 are the cell average values of the ALE-DG solution u_h and the quantities |K^n|, |K^n,1|, |K^n+1| are computed by (<ref>). We note that |K(t_n+1)|=|K^n+1|, since in two space dimensions, the TVD-RK2 method in Table <ref> solves the ODE (<ref>) exactly. 


The TVD-RK3 ALE-DG method:
For the TVD-RK3 method in Table <ref> the method (<ref>) becomes


    |K^n,1|=   |K^n|+ t(∇_ξ·[ A_K^n^-1(ω^n)])|K^n|, 
    |K^n,2|=   3/4|K^n|+1/4|K^n,1|+ t/4(∇_ξ·[ A_K^n+1^-1(ω^n+1)])|K^n,1|, 
    |K^n+1|=   1/3|K^n|+2/3|K^n,2|+2 t/3(∇_ξ·[ A_K^n+1/2^-1(ω^n+1/2)])|K^n,2|,


where t_n+1/2:=1/2(t_n+t_n+1) and K^n+1/2=K(t_n+1/2). In two and three space dimensions, the third order TVD-RK3 method in Table <ref> solves the ODE (<ref>) exactly and thus |K(t_n+1)|=|K^n+1|.
Next, we apply the test function v^*=1 in (<ref>) and obtain the scheme


    |K^n,1|u_K^n,1^n,1=   |K^n|u_K^n^n + t𝒢(u_h^n,*,1,J_K^n,t_n), 
    |K^n,2|u_K^n,2^n,2=   3/4|K^n|u_K^n^n+1/4|K^n,1|u_K^n,1^n,1+ t/4𝒢(u_h^n,1,*,1,J_K^n,1,t_n+1), 
    |K^n+1|u_K^n+1^n+1=   1/3|K^n|u_K^n^n+2/3|K^n,2|u_K^n,2^n,2+2 t/3𝒢(u_h^n,2,*,1,J_K^n,2,t_n+1/2),


where the values u_K^n^n, u_K^n,1^n,1, u_K^n,2^n,2 and u_K^n+1^n+1 are the cell average values of the ALE-DG solution u_h and the quantities |K^n|, |K^n,1|, |K^n,2|, |K^n+1| are computed by (<ref>).




  §.§.§ A quadrature rule to decompose the cell average value

In order to prove that the second and the third order fully-discrete ALE-DG methods satisfy the maximum principle, when the ALE-DG solution is revised by the limiter (<ref>), we proceed similar to the derivation by Zhang, Xia and Shu in <cit.>. For this task, we need to apply a special quadrature formula developed by Zhang, Xia and Shu to decompose the cell average values of the ALE-DG solution. In the following this quadrature formula is briefly reviewed.

First of all, it should be noted that in the implementation of the ALE-DG method the edge integrals are approximated  by a k+1-point 2k+1 accurate Gauss quadrature formula. Hence, we obtain

    ⟨g(ω,u_h^*,int_,u_h^*,ext_,J_K(t)(t)),1⟩ _∂
    ≈   ∑_β=1^k+1∑_ν=1^3σ_βg(ω_ν,β,u_ν,β^*,int_,u_ν,β^*,ext_,J_K(t)_F_^ν(t))ℓ_F_^ν,

where (t)= A_K(t)^-T_, _F_^ν(t)= A_K(t)^-T_F_^ν, F_^ν, ν =1,2,3, are the edges of the reference cell , the corresponding normals and lengths of the edges are _F_^ν  as well as  ℓ_F_^ν. Moreover, we denote by u_ν,β^*,int_ as well as u_ν,β^*,ext_ the values of the mapped ALE-DG solution u_h^*:=u_h∘χ_K(t) evaluated in the β-th Gauss quadrature point on the edge F_^ν. The corresponding quadrature weights for the interval [-1/2,1/2] are σ_β. Likewise, ω_ν,β are the values of the grid velocity ω evaluated in the β-th Gauss quadrature point on the edge F_^ν. We obtain for a constant c∈

    ⟨(ωc)· J_K(t)(t),1⟩ _∂=∑_β=1^k+1∑_ν=1^3σ_β(ω_ν,βc)·(J_K(t)_F_^ν(t))ℓ_F_^ν,

since the edge integrals are approximated by a 2k+1 accurate Gauss quadrature formula and the grid velocity belongs to the space  P^1(,^d) for all t∈[t_n,t_n+1] according to <Ref>. Thus, it follows for a constant c∈

    ∑_β=1^k+1∑_ν=1^3σ_βg(ω_ν,β,c,c,J_K(t)_F_^ν(t))ℓ_F_^ν
    
    =    ⟨(f(c)-ωc)·(J_K(t)(t)),1⟩ _∂=⟨g(ω,c)·(J_K(t)(t)),1⟩ _∂,

since the numerical flux has the property (P1). Therefore, the statement of <Ref> stays true when we approximate the edge integrals in (<ref>) by a 2k+1 accurate Gauss quadrature formula.

The cell average values need to be decomposed by a quadrature formula, which includes the Gauss quadrature points for the edges F_^ν, ν =1,2,3. A quadrature formula with this property has been developed by Zhang et al. in <cit.>. Let us assume that N is the smallest integer with 2N-3≥ k. The 3(N-1)(k+1)-point quadrature formula for triangular elements in <cit.>, has the properties:


  * The quadrature formula is exact for the integration of polynomials with degree k on a triangular element.

  * The quadrature points include the Gauss quadrature points for the edges F_^ν, ν =1,2,3.

  * All the quadrature weights are positive and the weights for the Gauss quadrature points are given by

    σ_βσ̂:=2/3σ_βσ̃,  σ̃=1/N(N-1),
      β=1,...,k+1,

where σ̃ corresponds to the 1-st and N-th Gauss-Lobatto quadrature weights for the interval [-1/2,1/2].

  * The quadrature formula has L:=3(N-2)(k+1) points in the interior of a triangular element.

This quadrature formula ensures that the cell average values can be written as follows

    u_K(t)=∑_ν=1^3∑_β=1^k+1σ_βσ̂u_ν,β^*,int_+∑_γ=1^Lσ_γu_γ^*,int_.

We denote by u_γ^*,int_, γ=1,...,L, the values of the mapped ALE-DG solution u_h^*:=u_h∘χ_K(t) evaluated in the quadrature points which are lying in the interior of the reference cell . The corresponding quadrature weights are denoted by σ_γ.






  §.§.§ The maximum principle for the TVD-RK2 ALE-DG method

In this section, we prove that the TVD-RK2 ALE-DG method satisfies the maximum principle, when the ALE-DG solution is revised by the bound preserving limiter (<ref>).

First of all, we show the property (L1). Therefore, similar to the process in <cit.>, we decompose the stages of the scheme (<ref>) in a sum of monotone increasing functions, which preserve constant states.
Therefore, it is convenient to use vector notations. In particular, we define the set

    𝕄:={ v=( v_1, v_2, v_3):  v_ν=(v_ν,1,...,v_ν,k+1)∈[m,M]^k+1, ∀ν=1,2,3}

with m and M given by (<ref>) and apply for any cell K(t)∈𝒯_(t) and all ν=1,2,3 the vector notations


    u_ν^*,int_:=(u_ν,1^*,int_,… u_ν,k+1^*,int_),    u_ν^*,ext_:=(u_ν,1^*,ext_,… u_ν,k+1^*,ext_),


    u^*,int_:=(u_1^*,int_,…,u_L^*,int_),


    u^*,int_:=( u_1^*,int_, u_2^*,int_, u_3^*,int_),    u^*,ext_:=( u_1^*,ext_, u_2^*,ext_, u_3^*,ext_).


Then, by applying the decomposition (<ref>) of the cell average values, the scheme  (<ref>) can be written as


    u_K^n,1^n,1=   ℒ(u^n,*,int_, u^n,*,int_, u^n,*,ext_,|K^n,1|,|K^n|,t_n), 
    u_K^n+1^n+1=   1/2ℋ(u^n,*,int_, u^n,*,int_,|K^n+1|,|K^n|) 
       +1/2ℒ(u^n,1,*,int_, u^n,1,*,int_, u^n,1,*,ext_,|K^n+1|,|K^n,1|,t_n+1),


where for all a∈[m,M]^L and b, c∈𝕄


    ℋ( a, b,|K_1|,|K_2|)
    :=   ∑_γ=1^Lσ_γ(1-1/|K_1|(|K_1|-|K_2|))a_γ
       +∑_β=1^k+1∑_ν=1^3σ_βσ̂(1-1/|K_1|(|K_1|-|K_2|))b_ν,β, 
       
    ℒ( a, b, c,|K_1|,|K_2|,t)
    :=   ∑_γ=1^Lσ_γ(1-1/|K_1|(|K_1|-|K_2|))a_γ
        +∑_β=1^k+1σ_βσ̂H_1(b_1,β,b_2,β,c_1,β,|K_1|,|K_2|,t)  
       +∑_β=1^k+1σ_βσ̂H_2(b_1,β,b_2,β,b_3,β,c_2,β,|K_1|,|K_2|,t) 
        +∑_β=1^k+1σ_βσ̂H_3(b_2,β,b_3,β,c_3,β,|K_1|,|K_2|,t)

with

    H_1(b_1,β,b_2,β,c_1,β,|K_1|,|K_2|,t)
    =   (1-1/|K_1|(|K_1|-|K_2|))b_1,β
        - t/σ̂|K_1|g(ω_1,β,b_1,β,c_1,β,J_K_2 A_K(t)^-T_F_^1)ℓ_F_^1
        + t/σ̂|K_1|g(ω_1,β,b_2,β,b_1,β,J_K_2 A_K(t)^-T_F_^1)ℓ_F_^1, 
       
    
    H_2(b_1,β,b_2,β,b_3,β,c_2,β,|K_1|,|K_2|,t)
    =   (1-1/|K_1|(|K_1|-|K_2|))b_2,β
        - t/σ̂|K_1|g(ω_1,β,b_2,β,b_1,β,J_K_2 A_K(t)^-T_F_^1)ℓ_F_^1
        - t/σ̂|K_1|g(ω_2,β,b_2,β,c_2,β,J_K_2 A_K(t)^-T_F_^2)ℓ_F_^2
       - t/σ̂|K_1|g(ω_3,β,b_2,β,b_3,β,J_K_2 A_K(t)^-T_F_^3)ℓ_F_^3, 
       
    
    H_3(b_2,β,b_3,β,c_3,β,|K_1|,|K_2|,t)
    =   (1-1/|K_1|(|K_1|-|K_2|))b_3,β
       - t/σ̂|K_1|g(ω_3,β,b_3,β,c_3,β,J_K_2 A_K(t)^-T_F_^3)ℓ_F_^3
       + t/σ̂|K_1|g(ω_3,β,b_2,β,b_3,β,J_K_2 A_K(t)^-T_F_^3)ℓ_F_^3.



In the Section <ref>, we mentioned that the result in <Ref> also holds when the edge integrals in (<ref>) are approximated by a 2k+1 accurate Gauss quadrature formula. Therefore for a vector ( c, c, c)∈[m,M]^L∪𝕄 with vector components given by the constant c∈[m,M], we obtain by the <Ref>

    ℒ( c, c, c,|K^n,1|,|K^n|,t_n)=c,


    1/2ℋ( c, c,|K^n+1|,|K^n|)+1/2ℒ( c, c, c,|K^n+1|,|K^n,1|,t_n+1)=c.

Henceforth, for the sake of simplicity, we will use the Lax-Friedrichs flux (<ref>) for the analysis. However, it should be noted that the techniques which are presented in this section can be also applied to any other monotone flux. The use of the Lax-Friedrichs flux ensures to prove the following lemmas.


Let |K^n|, |K^n,1|, |K^n+1| be given by (<ref>) and a∈[m,M]^L, b, c∈𝕄. Then, under the CFL constraint

    max{σ̂|(∇_ξ·[ A_K(t)^-1(ω(t))])||K(t)|+∑_ν=1^3λ^nℓ_F_^ν: t∈[t_n,t_n+1]} t/|K^n,1|≤σ̂

with λ^n given by (<ref>), it holds

    m≤ℒ( a, b, c,|K^n,1|,|K^n|,t_n)≤ M.




First of all, we apply the Lax-Friedrichs flux (<ref>) and rewrite the functions (<ref>), (<ref>) and (<ref>) as follows

    H_1(b_1,β,b_2,β,c_1,β,|K_1|,|K_2|,t) 
    
    =   (1-1/|K_1|(|K_1|-|K_2|)- t/σ̂|K_1|λ^nℓ_F_^1)b_1,β
       + t/σ̂|K_1|(g_+(ω_1,β,b_2,β,J_K_2_F_^1(t))+g_-(ω_1,β,c_1,β,J_K_2_F_^1(t)))ℓ_F_^1, 
       
        H_2(b_1,β,b_2,β,b_3,β,c_2,β,|K_1|,|K_2|,t)
    
    =   (1-1/|K_1|(|K_1|-|K_2|))b_2,β- t/σ̂|K_1|∑_ν=1^3g_+(ω_ν,β,b_2,β,J_K_2_F_^ν(t))ℓ_F_^ν
       + t/σ̂|K_1|g_-(ω_1,β,b_1,β,J_K_2_F_^1(t))ℓ_F_^1
    + t/σ̂|K_1|g_-(ω_2,β,c_2,β,J_K_2_F_^2(t))ℓ_F_^2
       + t/σ̂|K_1|g_-(ω_3,β,b_3,β,J_K_2_F_^3(t))ℓ_F_^3, 
       
        H_3(b_2,β,b_3,β,c_3,β,|K_1|,|K_2|,t) 
    
    =   (1-1/|K_1|(|K_1|-|K_2|)- t/σ̂|K_1|λ^nℓ_F_^3)b_3,β
       + t/σ̂|K_1|(g_+(ω_3,β,b_2,β,J_K_2_F_^3(t))+g_-(ω_3,β,c_3,β,J_K_2_F_^3(t)))ℓ_F_^3,

where _F_^ν(t):= A_K(t)^-T_F_^ν, ν=1,2,3.

Next, we observe that for all t∈[0,T], ν=1,2,3 and β=1,...,k+1

    |∂_b_2,βg_+(ω_ν,β,b_2,β,J_K(t)_F_^ν(t))|≤λ^n,

where λ^n is given by (<ref>). Hence, we obtain by (<ref>) as well as the CFL constraint (<ref>)
for  β=1,...,k+1 

    ∂_b_2,βH_2(b_1,β,b_2,β,b_3,β,c_2,β,|K^n,1|,|K^n|,t_n) 
    ≥     1- t/σ̂|K^n,1|(σ̂|∇_ξ·[ A_K^n^-1(ω^n)]||K^n|+∑_ν=1^3λ^nℓ_F_^ν)≥0,

and similar we get

    ∂_b_1,βH_1(b_1,β,b_2,β,c_1,β,|K^n,1|,|K^n|,t_n)≥0,


    ∂_b_3,βH_3(b_2,β,b_3,β,c_3,β,|K^n,1|,|K^n|,t_n)≥0.


In the following, we highlight by the symbolic notation ↑ that a function is increasing in the marked arguments. Likewise, we apply the notation ↑ to highlight that a function with vector arguments increases in each vector component. Then, it follows by (<ref>), (<ref>) and (<ref>)

    H_1(↑,↑,↑,|K^n,1|,|K^n|,t_n), H_2(↑,↑,↑,↑,|K^n,1|,|K^n|,t_n), H_3(↑,↑,↑,|K^n,1|,|K^n|,t_n),

since g_±(ω,↑,J_K(t)(t)) for all t∈[t_n,t_n+1] and all cells K(t)∈𝒯_(t).
Therefore, we obtain

    ℒ(↑,↑,↑,|K^n,1|,|K^n|,t_n),

since for all γ =1,...,L, it follows by (<ref>) and the CFL constraint (<ref>)

    ∂_a_γℒ( a, b, c,|K^n,1|,|K^n|,t_n)
    ≥σ_γ(1- t/|K^n,1||∇_ξ·[ A_K^n^-1(ω^n)]||K^n|)≥0.

Finally, the equation (<ref>) and (<ref>) provide the inequality (<ref>).



Let |K^n|, |K^n,1|, |K^n+1| be given by (<ref>) and a, a∈[m,M]^L, b,  b, c∈𝕄. Then, under the CFL constraint

    max{σ̂|∇_ξ·[ A_K(t)^-1(ω(t))]||K(t)|+∑_ν=1^3λ^nℓ_F_^ν: t∈[t_n,t_n+1]} t/|K^n+1|≤σ̂

with λ^n given by (<ref>), we have

    m≤1/2ℋ( a, b,|K^n+1|,|K^n|)+1/2ℒ( a, b, c,|K^n+1|,|K^n,1|,t_n+1)≤ M.



The equations (<ref>) and (<ref>) supply

    |K^n+1|-|K^n,1|= t/2(∇_ξ·[ A_K^n+1^-1(ω^n+1)]|K^n,1|-∇_ξ·[ A_K^n^-1(ω^n)]|K^n|).

Thus, the same procedure as in the proof of <ref> provides

    ℒ(↑,↑,↑,|K^n+1|,|K^n,1|,t_n+1)

by applying the identity (<ref>) and the CFL constraint (<ref>). Furthermore, by (<ref>) and (<ref>) follows

    |K^n+1|-|K^n|
    =
     t/2(∇_ξ·[ A_K^n^-1(ω^n)]|K^n|+∇_ξ·[ A_K^n+1^-1(ω^n+1)]|K^n,1|).

Hence, we obtain for all γ =1,...,L

    ∂_a_γℋ( a, b, |K^n+1|,|K^n|,t_n+1) 
    ≥   σ_γ(1- t/2|K^n+1|(∇_ξ·[ A_K^n^-1(ω^n)]|K^n|+|∇_ξ·[ A_K^n+1^-1(ω^n+1)]||K^n,1|))≥0,

by the identity (<ref>) and the CFL constraint (<ref>). In a similar way, it follows for all ν=1,2,3 as well as β =1,...,k+1

    ∂_b_ν,βℋ( a,b, |K^n+1|,|K^n|,t_n+1)≥0.

This ensures

    ℋ(↑,↑,|K^n+1|,|K^n|).

Therefore, the equation (<ref>), (<ref>) and (<ref>) supply the inequality (<ref>).


We note that the assumption (A3) for the mesh parameter provides for all t∈[0,T] and all cells K(t)∈𝒯_(t)

    |K(t)|≥πρ_K(t)^2≥π/κ^2h_K(t)^2≥π/τ^2κ^2h^2,

since ρ_K(t) denotes the radius of the largest ball contained in the cell K(t). Therefore, the CFL constraints (<ref>) and (<ref>) can be generalized as

    max{σ̂|∇_ξ·[ A_K(t)^-1(ω(t))]||K(t)|+∑_ν=1^3λ^nℓ_F_^ν: t∈[t_n,t_n+1]} t/h^2≤πσ̂/τ^2κ^2.

Now, we apply the generalized CFL constraint and the previous lemmas to prove the maximum principle for the second order fully-discrete ALE-DG method.


Suppose u^n,*,int_∈[m,M]^L u^n,*,int_, u^n,*,ext_∈𝕄. Furthermore, (A1) - (A3) and the CFL constraint (<ref>) are satisfied. Then, the solution u_h^n+1 of the second order fully-discrete ALE-DG method revised by Zhang, Xia and Shu's bound preserving limiter (<ref>) belongs to the interval [m,M].




Let K^n∈𝒯_(t_n) be an arbitrary cell. The main part of the proof follows in two steps.

Step 1.
Since, u^n,*,int_∈[m,M]^L u^n,*,int_, u^n,*,ext_∈𝕄 and the CFL constraint (<ref>) is satisfied, it follows u_K^n,1^n,1∈[m,M] by the inequality (<ref>) in <ref> and the equation (<ref>). Hence, the bound preserving limiter (<ref>) ensures u_h^n,1|_K^n,1∈[m,M].

Step 2.
In the first step, it has been shown that u_h^n,1|_K^n,1∈[m,M] when the bound preserving limiter (<ref>) was applied. Hence, u^n,1*,int_∈[m,M]^L, u^n,1,*,int_, u^n,1,*,ext_∈𝕄.

Thus, the inequality (<ref>) in <ref> and the equation (<ref>) supply u_K^n+1^n+1∈[m,M], since the CFL constraint (<ref>) is satisfied. Finally, the bound preserving limiter (<ref>) ensures that u_h^n+1|_K^n+1∈[m,M].

In a similar way, we proceed for any other cell in 𝒯_(t_n+1). Thus, it follows u_h^n+1∈[m,M].






  §.§.§ The maximum principle for the third order fully-discrete ALE-DG method

In this section, we show briefly how the result in <ref> can be extended to the third order TVD-RK3 ALE-DG method. We apply the decomposition (<ref>) of the cell average values and the vector notations (<ref>) to rewrite the scheme (<ref>) as

    u_K^n,1^n,1=   ℒ(u^n,*,int_, u^n,*,int_, u^n,*,ext_,|K^n,1|,|K^n|,t_n), 
    
    u_K^n,2^n,2=   3/4ℋ(u^n,*,int_, u^n,*,int_,|K^n,2|,|K^n|) 
       +1/4ℒ(u^n,1,*,int_, u^n,1,*,int_, u^n,1,*,ext_,|K^n,2|,|K^n,1|,t_n+1),
    
    u_K^n+1^n+1=   1/3ℋ(u^n,*,int_, u^n,*,int_,|K^n+1|,|K^n|)
         +2/3ℒ(u^n,2,*,int_, u^n,2,*,int_, u^n,2,*,ext_,|K^n+1|,|K^n,2|,t_n+1/2).


By the <Ref> we obtain the identities

    ℒ( c, c, c,|K^n,1|,|K^n|,t_n)=c,


    3/4ℋ( c, c,|K^n,2|,|K^n|)+1/4ℒ( c, c, c,|K^n,2|,|K^n,1|,t_n+1)=c,


    1/3ℋ( c, c,|K^n+1|,|K^n|)+2/3ℒ( c, c, c,|K^n+1|,|K^n,2|,t_n+1/2)=c,

where ( c, c, c)∈[m,M]^L∪𝕄 with vector components given by the constant c∈[m,M]. Moreover, like in the previous section, we apply the Lax-Friedrichs flux (<ref>). Then, by the same argumentation as in the proof of the  <ref> and <ref>, we obtain the following lemma.

Let |K^n|, |K^n,1|, |K^n,2|, |K^n+1| be given by (<ref>) and a, a∈[m,M]^L, b, c,b,c∈𝕄. Then, under the CFL constraint (<ref>), it holds

    m≤ℒ( a, b, c,|K^n,1|,|K^n|,t_n)≤ M,


    m≤3/4ℋ( a, b,|K^n,2|,|K^n|)+1/4ℒ( a, b, c,|K^n,2|,|K^n,1|,t_n+1)≤ M,


    m≤1/3ℋ( a, b,|K^n+1|,|K^n|)+2/3ℒ( a, b, c,|K^n+1|,|K^n,2|,t_n+1/2)≤ M.



This lemma provide the following analogue of the Theorem <ref> for the second order fully-discrete ALE-DG method. The proof follows similar to the proof of Theorem <ref>. Hence, it is skipped.


Suppose u^n,*,int_∈[m,M]^L, u^n,*,int_, u^n,*,ext_∈𝕄. Furthermore, (A1) - (A3) and the CFL constraint (<ref>) are satisfied. Then, the solution u_h^n+1 of the third order fully-discrete ALE-DG method revised by Zhang, Xia and Shu's bound preserving limiter (<ref>) belongs to the interval [m,M].


It is also possible to apply other TVD-RK methods like the five stage fourth order method of Spiteri and Ruuth <cit.> as time integrator in a fully-discrete ALE-DG method. Then, it can be proven by the techniques presented in this section that these fully-discrete ALE-DG methods satisfy the maximum principle.






§ NUMERICAL EXPERIMENTS


In this section, we demonstrate the performance of the ALE-DG method for conservation laws in two dimensions. In our simulation, the criss-triangular meshes are used. Furthermore, the third order TVD Runge-Kutta method is used in the first example (Example 5.1). In the other examples the five stage fourth order TVD Runge-Kutta method of Spiteri and Ruuth <cit.> is used for the time discretization. We observe that for this high order approximation in time our theoretical results hold numerically, too. In order to avoid complications with the stability of the explicit time integrator, we apply the suitable CFL condition dependent on equation (<ref>).

To verify our theoretical results, we present numerical simulations for a linear advection equation and Burgers' equation. Moreover, to highlight that the ALE-DG method can be also used for systems of conservation laws, we present a plain wave problem and a smooth vortex problem for the compressible Euler equations with a polytropic gas.

In all the numerical simulations, we consider two moving mesh scenarios. First a static uniform criss-triangular mesh with cell size h_0 is used. Next, a moving mesh with the grid point distribution

    x_j(t_n)= x_j(0) + 0.3sin(2π x_j(0)/x_r-x_l)sin(2π y_j(0)/y_r-y_l)sin(2π(t_n)/t_0),
       y_j(t_n)= y_j(0) + 0.2sin(2π x_j(0)/x_r-x_l)sin(2π y_j(0)/y_r-y_l)sin(4π(t_n)/t_0)

is used. In equation (<ref>) the points (x_j,y_j) are the vertices of the triangular mesh and t_0=√(10^2+5^2). The vertices (x_j(0),y_j(0)) at initial time are given by the same mesh as in the first moving mesh scenario. The grid point distribution (<ref>) has also been used by  Klingenberg <cit.> et al. and Persson et al. <cit.>. As an example, in Fig. <ref> we draw a typical mesh at t=0 and the deformed one at t=1.






  
Example 5.1  (Linear advection equation)
Here, we test the linear equation

    ∂_t u + ∂_x u+∂_y u=0,    (x,y)∈ [0,2]×[0,2],

with the periodic boundary condition and the initial condition is taken to be u_0(x,y)=1+0.5sin(π(x+y)). The exact solution is u(x,y,t)=u_0(x-t,y-t) at time t.

The numerical solution on the static uniform grid is u_h^S and the moving mesh solution is u_h^M. In Table <ref>, we show the L^2-errors and the rates of convergence of the numerical solutions u_h^S and u_h^M  for the advection equation at time t=1. In the computation, we used piecewise P^k polynomial spaces with k=1,2,3 on static and moving triangular meshes with cell size h_0. We observe that both u_h^S and u_h^M  have the optimal accuracy when P^k polynomial spaces with k=1,2,3 are applied, but they do not satisfy the discrete maximum principle. This was expected, since the bound preserving limiter in Section <ref> was not applied in the test case.

To verify the maximum principle, we run the same test case again, but this time with the bound preserving limiter. The results are given in Table <ref>, where u_h^M is the ALE-DG solution on the moving mesh with bound preserving limiter. We observe that the L^2-errors and the rates of convergence for the moving mesh solution u_h^M are not affected by the use of the bound preserving limiter and thus the method still has optimal accuracy. Furthermore, we  can see that the numerical solution u_h^M is limited in the same range [0.5,1.5] as the initial data.










  
Example 5.2  (Burgers' equation)
Next, we investigate the Burgers' equation

    ∂_t u + ∂_x (u^2/2)+∂_y (u^2/2)=0,    (x,y)∈ [0,2]×[0,2].

In our simulation, the periodic boundary condition is used and the initial condition  is also taken to be u_0(x,y)=1+0.5sin(π(x+y)). We compute this example up to time t=0.1 before the shock front has been developed in the numerical solution.

In Table <ref>, the L^2-errors and the rates of convergence for the numerical solutions u_h^S and u_h^M, ũ_h^M are presented. These functions are computed by the ALE-DG method with P^k, k=1,2,3, polynomial spaces. As in the previous example u^S_h and u^M_h are the numerical solutions of the ALE-DG method on the static uniform mesh and on the moving mesh. The solution ũ_h^M is the moving mesh ALE-DG solution revised by the bound preserving limiter. We observe that the optimal accuracy is obtained for u_h^M,u_h^M and u_h^S in both moving mesh scenarios. Furthermore, maximum and minimum values of u_h^M are limited in the same range [0.5,1.5]  as the initial data when the bound preserving limiter is applied in the ALE-DG method on the moving grid.


To show that our proposed ALE-DG methods can handle the problem with shocks, we show the numerical solutions u_h^S and u_h^M of Burgers' equation at time t=0.45 with piecewise P^1 polynomial approximation in Fig. <ref>. Here, we use the slope limiter developed by Cockburn et al. in <cit.>. From the results, it can be seen that the ALE-DG methods can capture the shocks well for the Burgers' equation on both static and moving meshes.






  
Example 5.3  (Compressible Euler equations)
We consider the two dimensional compressible Euler equations of gas dynamics for a polytropic gas

    ∂_t𝐔+∇·𝐅(𝐔)=0,    (x,y)∈ [x_l,x_r]×[y_l,y_r]⊂^2,

with

    𝐔=(ρ,ρ u,ρ v,E)^T,  𝐅(𝐔)=[ρ u,ρ𝐮⊗𝐮+p I,(E+p)𝐮]^T.

Here, ρ is the density, 𝐮=(u,v)^T is the velocity field and E is the total energy. Moreover, the adiabatic constant of air  γ=1.4 is used, the pressure is given by p=(γ-1)(E-1/2ρ u^2) and I is the identity matrix.
In our simulation, we test a plain wave problem and a smooth vortex problem.

First, we consider the plain wave problem and choose the domain related parameter in (<ref>) as x_l=y_l=0 and x_r=y_r=2. The problem has the initial data

    (ρ, u, v, p)^T=(1+0.5sin(π(x+y)),1,1,1)^T

and is investigated with the periodic boundary condition. The results in Table <ref> show the L^2-errors and the rates of convergence of the density ρ_h given by the ALE-DG method with P^k, k=1,2,3, polynomial spaces. The numerical solutions of the ALE-DG method on the static uniform mesh and the moving mesh are ρ^S_h and ρ^M_h. The numerical results show that we can obtain the optimal accuracy on both meshes.



Next, we consider the smooth vortex problem and choose the domain related parameter in (<ref>) as x_l=y_l=0, x_r=20 and y_r=15. This problem was also presented by Persson et al. <cit.> and the initial condition is

    ρ= ρ_0(1-α e^r)^1/γ-1,    p= p_0(1-α e^r)^γ/γ-1,
        u=(u,v)^T=(u_0cos(θ),v_0sin(θ))^T+ϵ/2π r_0e^0.5r(-u_0(y-y_0),v_0(x-x_0))^T,

where (ρ_0,u_0,v_0,p_0)^T=(1,1,1,1)^T, θ=arctan(0.5), (x_0,y_0)=(5,5), ϵ=0.3, r_0=1.5, r=(1-(x-x_0)^2-(y-y_0)^2)/r_0^2 and α=(γ-1)ϵ^2/8γπ^2. We test the problem up to time t=√(10^2+5^2) with the Direchlet boundary condition. The ALE-DG method with P^k,k=1,2,3 approxiation is used  to solve the problem  on static uniform triangular meshes with cell size h_0 and on moving meshes with the grid point distribution (<ref>) as before. The L^2-errors and the rates of convergence for the numerical solutions of the density ρ_h^S,ρ_h^M and the pressure p_h^M, p_h^S are shown in Table <ref>. We see that the numerical solutions are optimally accurate in both moving mesh scenarios.







  
Example 5.4  (Constant state preservation)
The previous examples show that the ALE-DG method on moving meshes maintains the high order accuracy as the DG method on static meshes. The ability of the ALE-DG method to preserve constant states needs to be investigated, too. For this reason the linear advection equation (<ref>)  and the Burgers' equation (<ref>) are considered with the constant initial condition u_0=1. We solve these initial value problems with the ALE-DG method on moving triangular meshes with the grid point distribution (<ref>). In Table <ref> the results of the computations are listed and it can be seen that the ALE-DG method numerically satisfies the GCL. This result was expected, since we used a time discretization with an order greater than two and in the Section <ref> it has been proven that a time discretization of this type is enough to ensure that the method preserves constant states.




To further show the D-GCL of ALE-DG methods, we adopted the meshes 𝒯_1 and 𝒯_2 of 32 cells in  Fig. <ref> with recursive refinement as the initial and final meshes. For the linear advection equation with the constant initial condition u_0(x) = 1, we show the D-GCL errors at time T=1.0 in Table <ref> by forward Euler, TVD-RK2 and TVD-RK3 methods respectively. Here, P^1 piecewise polynomial space is used in the ALE-DG method. We take time step size t=h_0/max(|ω|) with ω=((𝐱_1-𝐱_2)/T,(𝐲_1-𝐲_2)/T) and (𝐱_1,𝐲_1),(𝐱_2,𝐲_2) are vertices of meshes 𝒯_1 and 𝒯_2. The numerical results are consistent with the analysis on the D-GCL of ALE-DG methods.




























§ CONCLUSIONS

In this paper, an ALE-DG method to solve conservation laws in several space dimensions on moving simplex meshes has been developed and analyzed. We began the paper with an analysis of the semi-discrete ALE-DG method and proved the L^2-stability. Moreover, we presented a suboptimal a priori error estimate with respect to the L^∞(0,T;L^2(Ω))-norm, where the suboptimality refers to the approximation properties of the discrete space.

Afterward, the fully-discrete ALE-DG method was investigated. In the context of Total-variation-diminishing Runge-Kutta methods, a relationship between the spatial dimension and the discrete geometric conservation law was elaborated. 
Furthermore, in two dimensions, second and third order fully-discrete ALE-DG methods were presented. We proved that these methods satisfy the maximum principle when the bound preserving limiter developed by Zhang, Xia and Shu <cit.> is applied. In a future work, it would be worthwhile to investigate if these methods are positive preserving when they are applied to the compressible Euler equations.

Beside our theoretical investigations, several numerical test examples for two moving mesh scenarios have been presented. These examples support our theoretical results and show that the ALE-DG method is numerically stable and uniformly high order accurate. In particular, the two test examples for the compressible Euler equations support the expectation that the ALE-DG method can be also applied to systems of conservation laws even when the development and analysis in this paper has been focused on scalar conservation laws in several space dimensions.

It should be mentioned that we did not use a moving mesh methodology in the numerical examples. The grid point distribution was specified for the calculations. The development of a suitable moving mesh methodology for our ALE-DG method is also a project for a future work.




99

Bellman R. Bellman. Introduction to matrix analysis, volume 19 of Classics in Applied Mathematics. SIAM, Philadelphia, PA, 2nd Edition, 1987.

Boscheri2017 W. Boscheri and M. Dumbser. Arbitrary-Lagrangian-Eulerian discontinuous Galerkin schemes with a posteriori subcell finite volume limiting on moving unstructured meshes. J. Comput. Phys. 346 (2017), 449-479.

CiarletNF P. G. Ciarlet. Linear and nonlinear functional analysis with applications, Volume 130, SIAM, Philadelphia, PA, 2013.

CiarletFE P. G. Ciarlet. The finite element method for elliptic problems, Volume 40 of Classics in Applied Mathematics, SIAM, Philadelphia, PA, 2002.

Cockburn1998 B. Cockburn and C.-W. Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws V, J. Comput. Phys. 141 (1998), 199-224.

Cockburn2001 B. Cockburn and C.-W. Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (2001), 173-261.

Donea2004 J. Donea, A. Huerta, J. P. Ponthot and A. Rodríguez-Ferran. In Encyclopedia of Computational Mechanics, E. Stein, R. De Borst and T. J.R. Hughes (Eds.), Volume. 1: Fundamentals., Chapter 14: Arbitrary Lagrangian-Eulerian Methods, Wiley, 2004.

Etienne2009 S. Étienne, A. Garon and D. Pelletier. Perspective on the geometric conservation law and finite element methods for ALE simulations of incompressible flow, J. Comput. Phys. 228 (2009), 2313-2333.

Farhat2001 C. Farhat, P. Geuzaine and C. Grandmont. The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids, J. Comput. Phys. 174 (2001), 669-694.

Gottlieb1998 S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta schemes, Math. Comp. 67 (1998), 73-85.

Gottlieb2001 S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods. SIAM review, 43 (2001), 89-112.

Farhat2000 H. Guillard and C. Farhat. On the significance of the geometric conservation law for flow computations on moving meshes, Comput. Method. Appl. M. 190 (2000), 1467-1482.



Jiang1994MC G.-S. Jiang and C.-W. Shu. On a cell entropy inequality for discontinuous Galerkin methods, Math. Comp. 62 (1994), 531-538.

Klingenberg2016 C. Klingenberg, G. Schnücke and Y. Xia. Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for conservation laws: analysis and application in one dimension, Math. Comp. 86 (2017), 1203-1232.

KlingenbergHJ2016HJ C. Klingenberg, G. Schnücke and Y. Xia. An arbitrary Lagrangian-Eulerian local discontinuous Galerkin method for Hamilton-Jacobi equations, J. Sci. Comput., 73 (2017), 906-942.

Kopriva2006 D. A. Kopriva. Metric identities and the discontinuous spectral element method on curvilinear meshes. J. Sci. Comput., 26 (2006), 301-327.

Kopriva2016 D. A. Kopriva, A. R. Winters, M. Bohm and G. J. Gassner. A provably stable discontinuous Galerkin spectral element approximation for moving hexahedral meshes. Comput. Fluids 139 (2016), 148-160.

Kruzkov1970 S. N. Kružkov. First order quasilinear equations in several independent variables, Math. USSR-Sbornik, 10 (1970), 217-243.

Kuzmin2010 D. Kuzmin. A guide to numerical methods for transport equations, Lecture, Friedrich-Alexander-Universität Erlangen-Nürnberg (2010).

Kuzmin2001 D. Kuzmin and S. Turek. Flux correction tools for finite elements, J. Comput. Phys. 175 (2002), 525-558.

Lesoinne1996  M. Lesoinne and C. Farhat. Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aeroelastic computations, Comput. Method. Appl. M. 134 (1996), 71-90.

Lombard1979 C. K. Lombard and P. D. Thomas. Geometric conservation law and its application to flow computations on moving grids, AIAA J., 17 (1979), 1030-1037.

Lomtev1999JCP I. Lomtev, R. M. Kirby and G. E. Karniadakis. A discontinuous Galerkin ALE method for compressible viscous flows in moving domains, J. Comput. Phys. 155 (1999), 128-159.

Luo2015 J. Luo, C.-W. Shu, and Q. Zhang. A priori error estimates to smooth solutions of the third order Runge –Kutta discontinuous Galerkin method for symmetrizable systems of conservation laws, ESAIM: M2AN 49 (2015), 991-1018.


Mavriplis2006 D. J. Mavriplis and Z. Yang. Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes, J. Comput. Phys. 213 (2006), 557-573.

Minoli2011
C. A. Acosta Minoli and D. A. Kopriva. Discontinuous Galerkin spectral element approximations on moving meshes,  J. Comput. Phys. 230 (2011), 1876-1902.

Nguyen2010JFS V. T. Nguyen. An arbitrary Lagrangian-Eulerian discontinuous Galerkin method for simulations of flows over variable geometries, J. Fluid. Struct. 26 (2010), 312-329.

Osher1985 S. Osher. Convergence of generalized MUSCL schemes, SIAM J. Numer. Anal. 22 (1985), 947-961.

Persson2009 P. O. Persson, J. Bonet, J. Peraire. Discontinuous Galerkin solution of the Navier-Stokes equations on deformable domains, Comput. Method. Appl. M. 198 (2009), 1585-1595.

Reed1973 W.H. Reed, T.R. Hill. Triangular mesh method for the neutron transport equation, Technical report LA-UR-73-479, Los Alamos Scientific Laboratory, Los Alamos, NM, 1973.

Robinson1991 B. A. Robinson, H. T. Yang and J. T. Batina. Aeroelastic analysis of wings using the Euler equations with a deforming mesh. J. Aircraft 28 (1991), 781-788.

Shu C.-W. Shu. Total-variation-diminishing time discretizations, SIAM J. Sci. Stat. Comput. 9 (1988), 1073-1084.

Spiteri R. J. Spiteri and S. J. Ruuth. A new class of optimal high-order strong-stability-preserving time discretization methods, SIAM J. Numer. Anal. 40 (2002), 469-491.

Persson2015  L. Wang and P. O. Persson. High-order discontinuous Galerkin simulations on moving Domains using ALE formulations and local remeshing and projections, 53rd AIAA Aerospace Sciences Meeting, AIAA SciTech Forum (AIAA 2015-0820).
Available via <http://dx.doi.org/10.2514/6.2015-0820 >

Xu2007
Y. Xu and C.-W. Shu. Error estimates of the semi-discrete local discontinuous Galerkin method for nonlinear convection–diffusion and KdV equations, Comput. Methods Appl. Mech. Eng., 196 (2007), 3805-3822.

Xu2014
Z. Xu. Parametrized maximum principle preserving flux limiters for high order schemes solving hyperbolic conservation laws: one-dimensional scalar problem, Math. Comp. 83 (2014), 2213-2238.

QZhang2004
Q. Zhang and C.-W. Shu. Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws, SIAM J. Numer. Anal. 42 (2004), 641-666.

QZhang2006
Q. Zhang and C.-W. Shu. Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin method for symmetrizable systems of conservation law, SIAM J. Numer. Anal. 44 (2006), 1703-1720.

QZhang2010
Q. Zhang and C.-W. Shu. Stability analysis and a priori error estimates of the third order explicit Runge–Kutta discontinuous Galerkin method for scalar conservation laws, SIAM J. Numer. Anal., 48 (2010), 1038-1063.

QZhang2014
Q. Zhang and Chi-Wang Shu. Error estimates for the third order explicit Runge-Kutta discontinuous Galerkin method for a linear hyperbolic equation in one-dimension with discontinuous initial data, Numer. Math. 126 (2014), 703-740.

Zhang2010 X. Zhang and C.-W. Shu. On maximum principle-satisfying high order schemes for scalar conservation laws, J. Comput. Phys. 229 (2010), 3091-3120.

Zhang2012 X. Zhang, Y. Xia and C.-W. Shu. Maximum-principle-satisfying and positivity-preserving high order discontinuous Galerkin schemes for conservation laws on triangular meshes, J. Sci. Comput. 50 (2012), 29-62.

Zhou2018 L. Zhou, Y. Xia and C-W. Shu. Stability analysis and error estimates of arbitrary Lagrangian-Eulerian discontinuous Galerkin method coupled with Runge-Kutta time-marching for linear conservation laws, ESAIM: Mathematical Modelling and Numerical Analysis, to appear.



