1 Introduction
Finite volume schemes based on exact or approximate Riemann solvers are used for solving hyperbolic conservation laws like the Euler equations governing compressible flows. These schemes are able to compute discontinuous solutions in a stable manner since they have implicit dissipation built into them due to the upwind nature of the schemes. Higher order schemes are constructed following a reconstruction approach combined with a high order time integration scheme. Discontinuous Galerkin methods can be considered as higher order generalizations of finite volume methods which also make use of Riemann solver technology but do not need a reconstruction step since they evolve a polynomial solution inside each cell. While these methods are formally high order accurate on smooth solutions, they can still introduce too much numerical dissipation in some situations. Springel springel2010pur gives the example of a KelvinHelmholtz instability in which adding a large constant velocity to both states leads to suppression of the instability due to excessive numerical dissipation. This behaviour is attributed to the fact that fixed grid methods based on upwind schemes are not Galilean invariant. Upwind schemes, even when they are formally high order accurate, are found to be too dissipative when applied to turbulent flows Johnsen:2010:AHM:1670833.1670927 since the numerical viscosity can overwhelm the physical viscosity.
For the linear convection equation
, the first order upwind scheme has the modified partial differential equation
which shows that the numerical dissipation is proportional to
which is the wave speed. In case of Euler equations simulated with a Riemann solver, e.g., the Roe scheme, the wave speeds are related to the eigenvalues of the flux Jacobian and the numerical dissipation would be proportional to the absolute values of the eigenvalues, e.g.,
where is the fluid velocity and is the sound speed. This type of numerical viscosity is not Galilean invariant since the fluid velocity depends on the coordinate frame adopted for the description of the flow. Adding a large translational velocity to the coordinate frame will increase the numerical viscosity and reduce the accuracy of the numerical solution. Such high numerical viscosity can be eliminated or minimized if the grid moves along with the flow as in Lagrangian methods doi:10.1137/0731002 ; Carre:2009:CLH:1552584.1553086 ; Maire20092391 . However pure Lagrangian methods encounter the issue of large grid deformations that occur in highly sheared flows as in the KelvinHelmholtz problem requiring some form of remeshing. A related approach is to use arbitrary LagrangianEulerian approach HIRT1974227 ; doneaale where the mesh velocity can be chosen to be close to the local fluid velocity but may be regularized to maintain the mesh quality. Even in the ALE approach, it may be necessary to perform some local remeshing to prevent the grid quality from degrading. In springel2010pur , the mesh is regenerated after every time step based on a Delaunay triangulation, which allows it to maintain good mesh quality even when the fluid undergoes large shear deformation. However these methods have been restricted to second order accuracy as they rely on unstructured finite volume schemes on general polygonal/polyhedral cells, where achieving higher order accuracy is much more difficult compared to structured grids.Traditionally, ALE methods have been used for problems involving moving boundaries as in wing flutter, store separation and other problems involving fluid structure interaction LOMTEV1999128 ,Venkatasubban19951743 ,wang2015aa , Boscheri2017 , Gaburro2018 . In these applications, the main reason to use ALE is not to minimize the dissipation in upwind schemes but to account for the moving boundaries, and hence the grid velocities are chosen based on boundary motion and with a view to maintain good mesh quality. Another class of methods solve the PDE on moving meshes where the mesh motion is determined based on a monitor function which is designed to detect regions of large gradients in the solution, see ammmbook and the references therein. These methods achieve automatic clustering of grid points in regions of large gradients. ALE schemes have been used to compute multimaterial flows as in Luo2004304
, since they are useful to accurately track the material interface. The mesh velocity was chosen to be equal to the contact speed but away from the material contact, the velocity was chosen by linear interpolation and was not close to Lagrangian. There are other methods for choosing the mesh velocity which have been studied in
Boscheri2013 ; Lohner2004 . LaxWendroff type ALE schemes for compressible flows have been developed in Liu20098872 . Finite volume schemes based on ADER approach have been developed on unstructured grids Boscheri_Dumbser_2013 , Boscheri201648 , boscheriIJNMF2014 . The theoretical analysis of ALEDG schemes in the framework of RungeKutta time stepping for conservation laws has been done in gero .In the present work, we consider only the one dimensional problem in order to set down the fundamental principles with which in an upcoming work, we shall solve the multidimensional problem. The numerical method developed here will be usable in the multiple dimensions, but additional work is required in multiple dimensions to maintain a good mesh quality under fluid flow deformation. We develop an explicit discontinuous Galerkin scheme that is conservative on moving meshes and automatically satisfies the geometric conservation law. The scheme is a single step method which is achieved by using a predictor computed from a RungeKutta scheme that is local to each cell in the sense that it does not require any data from neighbouring cells and belongs to the class of schemes called ADER method. Due to the single step nature of the scheme, the TVD limiter has to be applied only once in each time step unlike in multistage RungeKutta schemes where the limiter is applied after each stage update. This nature of the ADER scheme can reduce its computational expense especially in multidimensional problems and while performing parallel computations. The mesh velocity is specified at each cell face as the local velocity with some smoothing. We analyze the positivity of the first order scheme using Rusanov flux and derive a CFL condition. The scheme is shown to be exact for steady moving contact waves and the solutions are invariant to the motion of the coordinate frame. Due to Lagrangian nature, the Roe scheme does not require any entropy fix. However, we identify the possibility of spurious contact waves arising in some situations. This is due to the vanishing of the eigenvalue corresponding to the contact wave. While the cell averages are well predicted, the higher moments of the solution can be inaccurate. This behaviour of Lagrangian DG schemes does not seem to have been reported in the literature. We propose a fix for the eigenvalue in the spirit of the entropy fix of Harten
Harten1983 that prevents the spurious contact waves from occurring in the solution. The methodology developed here will be extended to multidimensional flows in a future work with a view towards handling complex sheared flows.The rest of the paper is organized as follows. Section (2) introduces the Euler equation model that is used in the rest of the paper. In section (3), we explain the derivation of the scheme on a moving mesh together with the quadrature approximations and computation of mesh velocity. The computation of the predicted solution is detailed in section (4). The TVD type limiter is presented in section (5) for a nonuniform mesh, section (6) shows the positivity of the first order scheme and section (7) shows the preservation of constant states. The grid coarsening and refinement strategy is explained in section (8) while section (9) presents a series of numerical results.
2 Euler equations
The Euler equations model the conservation of mass, momentum and energy, and can be written as a system of coupled partial differential equations laws of the form
(1) 
where
is called the vector of
conserved variables and are the corresponding fluxes given byIn the above expressions, is the density, is the fluid velocity, is the pressure and is the total energy per unit volume, which for an ideal gas is given by , with being the ratio of specific heats at constant pressure and volume, and is the enthalpy. The Euler equations form a hyperbolic system; the flux Jacobian
has real eigenvalues and linearly independent eigenvectors. The eigenvalues are
where is the speed of sound and the corresponding right eigenvectors are given by(2) 
The hyperbolic property implies that can be diagonalized as where is the matrix formed by the right eigenvectors as the columns and is the diagonal matrix of eigenvalues.
3 Discontinuous Galerkin method
3.1 Mesh and solution space
Consider a partition of the domain into disjoint cells with the ’th cell being denoted by . As the notation shows, the cell boundaries are time dependent which means that the cell is moving in some specified manner. The time levels are denoted by with the time step . The boundaries of the cells move with a constant velocity in the time interval given by
which defines a cell in spacetime as shown in figure (1). The algorithm to choose the mesh velocity is explained in a later section. The location of the cell boundaries is given by
Let and denote the center of the cell and its length, i.e.,
Let be the continuous linear interpolation of the mesh velocity which is given by
We will approximate the solution of the conservation law by piecewise polynomials which are allowed to be discontinuous across the cell boundaries as shown in figure (2). For a given degree , the solution in the ’th cell is given by
where are the degrees of freedom associated with the ’th cell. The basis functions are defined in terms of Legendre polynomials
where is the Legendre polynomial of degree . The above definition of the basis functions implies the following orthogonality property
(3) 
We will sometimes also write the solution in the ’th cell in terms of the reference coordinates as
and we will use the same notation to denote both functions.
3.2 Derivation of the scheme
In order to derive the DG scheme on a moving mesh, let us introduce the change of variable given by
(4) 
For any , we now calculate the rate of change of the ’th moment of the solution starting from
wherein we used the change of variables given by (4). But we also have the inverse transform
and hence
Using the above relations, we can easily show that
Moreover
since is linear in and hence is constant inside each cell. Hence the ’th moment evolves according to
where we have transformed back to the physical coordinates and made use of the conservation law (1) to replace the time derivative of the solution with the flux derivative. Define the ALE flux
(5) 
Performing an integration by parts in the variable, we obtain
where we have introduced the numerical flux
which provides an approximation to the ALE flux, see Appendix. Integrating over the time interval and using (3), we obtain
The above scheme has an implicit nature since the unknown solution appears on the right hand side integrals whereas we only know the solution at time . In order to obtain an explicit scheme, we assume that we have available with us a predicted solution in the time interval , which is used in the time integrals to obtain an explicit scheme. Moreover, the integrals are computed using quadrature in space and time leading to the fully discrete scheme
where are weights for time quadrature and are weights for spatial quadrature. For the spatial integral, we will use point Gauss quadrature. For the time integral we will use midpoint rule for and two point Gauss quadrature for . Since the mesh is moving, the spatial quadrature points depend on the quadrature time though this is not clear from the notation. In practice, the integrals are computed by mapping the cell to the reference cell, and the basis functions and its derivatives are also evaluated on the reference cell. The quadrature points in the reference cell are independent of time due to the linear mesh evolution.
3.3 Mesh velocity
The mesh velocity must be close to the local fluid velocity in order to have a Lagrangian character to the scheme. Since the solution is discontinuous, there is no unique fluid velocity at the mesh boundaries. Some researchers, especially in the context of Lagrangian methods, solve a Riemann problem at the cell face to determine the face velocity. Since we use an ALE formulation, we do not require the exact fluid velocity which is anyway not available to use since we only have a predicted solution. Following the exact trajectory of the fluid would also lead to curved trajectories for the grid point, which is an unnecessary complication. In our work, we make two different choices for the mesh velocities:

The first choice is to take an average of the two velocities at every face In the numerical results, we refer to this as ADG.

The second choice is to solve a linearized Riemann problem at the face at time . In the numerical results, we refer to this as RDG. For simplicity of notation, let the solution to the left of the face be represented as and the solution to the right be represented as , then
We will also perform some smoothing of the mesh velocity, e.g., the actual face velocity is computed from
Note that our algorithm to choose the mesh velocity is very local and hence easy and efficient to implement as it does not require the solution of any global problems. In Springel springel2010pur , the mesh velocity is adjusted so that the cells remain nearly isotropic which leads to smoothly varying cell sizes. Such an approach leads to many parameters that need to be selected and we did not find a good way to make this choice that works well for a range of problems. Instead, we will make use of mesh refinement and coarsening to maintain the quality of cells, i.e., to prevent very small or large cells from occurring in the grid. The use of a DG scheme makes it easy to perform such local mesh adaptation without loss of accuracy.
Remark 1
Consider the application of the proposed ALEDG scheme to the linear advection equation . In this case the mesh velocity is equal to the advection velocity , i.e., the cells move along the characteristics. This implies that the ALE flux and also the numerical flux . Thus the DG scheme reduces to
so that the solution at time has been advected exactly to the solution at time . Note that there is no time step restriction involved in this case and the accuracy of the predicted solution is also not relevant. If the initial condition has a discontinuity coinciding with a cell face, then the scheme advects the discontinuity exactly without any diffusion.
4 Computing the predictor
The predicted solution is used to approximate the flux integrals over the time interval and the method to compute this must be local, i.e., it must not require solution from neighbouring cells. Several methods for computing the predictor have been reviewed in Gassner20114232 . The simplest approach is to use a Taylor expansion in space and time. Since the cells are moving, the Taylor expansion has to be performed along the trajectory of the mesh motion. For a second order scheme, an expansion retaining only linear terms in and is sufficient. Consider a quadrature point ; the Taylor expansion of the solution around the cell center and time level yields
and the predicted solution is given by truncating the Taylor expansion at linear terms, leading to
Using the conservation law, the time derivative is written as so that the predictor is given by
(7) 
The above predictor is used for the case of polynomial degree . This procedure can be extended to higher orders by including more terms in the Taylor expansion but the algebra becomes complicated. Instead we will adopt the approach of continuous explicit RungeKutta (CERK) schemes Owren:1992:DEC:141072.141096 to approximate the predictor.
Let us choose a set of distinct nodes, e.g., GaussLegendre or GaussLobatto nodes, which uniquely define the polynomial of degree . These nodes are moving with velocity , so that the time evolution of the solution at node is governed by
wherein we have made use of the PDE to write the time derivative in terms of spatial derivative of the flux. This equation is solved with initial condition
Using a RungeKutta scheme of sufficient order (see Appendix), we will approximate the solution at these nodes as
where , is the stage time and are certain polynomials related to the CERK scheme and given in the Appendix. Note that we are evolving the nodal values but the computation of requires the modal representation of the solution in order to calculate spatial derivative of the solution.
Once the predictor is computed as above, it must be evaluated at the quadrature point as follows. For each time quadrature point ,

Compute nodal values ,

For each , convert the nodal values to modal coefficients ,

Evaluate predictor
The conversion from nodal to modal values is accomplished through a Vandermonde matrix of size which is the same for every cell and can be inverted once before the iterations start. The predictor is also computed at the cell boundaries using the above procedure. Figure (3) and (4) show the quadrature points used in the second, third and fourth order scheme. For the second order scheme, the values at and points are obtained from the predictor based on Taylor expansion as given in equation (7). For third and fourth order schemes, the nodal values are evolved forward in time by the CERK scheme and evaluated at the points. The point values are converted to modal coefficients using which the solution at the points is computed.
Remark 2
By looping over each cell in the mesh, the predicted solution is computed in each cell and the cell integral in (3.2) is evaluated. The trace values and at the points needed for quadrature in time are computed and stored. These are later used in a loop over the cell faces where the numerical flux is evaluated. Thus the algorithm is easily parallelizable on multiple core machines and/or using threads.
5 Control of Oscillation by Limiting
High order schemes for hyperbolic equations suffer from spurious numerical oscillations when discontinuities or large gradients are present in the solution which cannot be accurately resolved on the mesh. In the case of scalar problems, this is a manifestation of loss of TVD property and hence limiters are used to satisfy some form of TVD condition. In the case of DG schemes, the limiter is used as a post processor which is applied on the solution after the time update has been performed. If the limiter detects that the solution is oscillatory, then the solution polynomial is reduced to at most a linear polynomial with a limited slope. In the present scheme, the limiter is applied after the solution is updated from time to time , i.e., the solution obtained from (3.2) is postprocessed by the limiter. Since the mesh is inherently nonuniform due to it being moved with the flow, we modify the standard TVD limiter to account for this nonuniformity. Also, since we are solving a system of conservation laws, the limiter is applied on the local characteristic variables which gives better results than applying it directly on the conserved variables Cockburn:1989:TRL:69978.69982 .
The solution in cell can be written as ^{1}^{1}1We suppress the time variable for clarity of notation.
where is the cell average value and is proportional to the derivative of the solution at the cell center. Let , denote the matrix of right and left eigenvectors evaluated at the cell average value which satisfy , and the right eigenvectors are given in equation (2). The local characteristic variables are defined by and . We first compute the limited slope of the characteristic variables from
where the minmod function is defined by
If then we retain the solution as it is, and otherwise, the solution is modified to
This corresponds to a TVD limiter which is known to lose accuracy at smooth extrema cockburn1989 since the minmod function returns zero slope at local extrema. The TVB limiter corresponds to replacing the minmod limiter function with the following function
(8) 
where the parameter
is an estimate of the second derivative of the solution at smooth extrema
cockburn1989 , and has to be chosen by the user.6 Positivity property
The solutions of Euler equations are well defined only if the density and pressure are positive quantities. This is not a priori guaranteed by the DG scheme even when the TVD limiter is applied. In the case of RungeKutta DG schemes, a positivity limiter has been developed in Zhang:2010:PHO:1864819.1865066 which preserves accuracy in smooth regions. This scheme is built on a positive first order finite volume scheme. Consider the first order version of the ALEDG scheme which is a finite volume scheme given by
(9) 
The only degree of freedom is the cell average value and the solution is piecewise constant. We will analyze the positivity of this scheme for the case of Rusanov flux which is given in A.1. The update equation can be rewritten as
From the definition of the Rusanov flux formula, we can easily see that^{2}^{2}2We drop the superscript in some of these expressions.
Consider the first component of
Consider the first component of
The pressure corresponding to is
and similarly the pressure corresponding to is nonnegative. Hence if the coefficient in term is positive, then the scheme is positive. This requires the CFL condition
The time step will also be restricted to ensure that the cell size does not change too much in one time step. If we demand that the cell size does not change by more than a fraction , then we need to ensure that the time step satisfies
Combining the previous two conditions, we obtain the following condition on the time step
(10) 
We can now state the following result on the positivity of the first order finite volume scheme on moving meshes.
Theorem 6.1
Remark 3
In this work we have not attempted to prove the positivity of the scheme for other numerical fluxes. We also do not have a proof of positivity for higher order version of the scheme. In the computations, we use the positivity preserving limiter of Zhang:2010:PHO:1864819.1865066 which leads to robust schemes which preserve the positivity of the cell average value in all the test cases.
7 Preservation of constant states
An important property of schemes on moving meshes is their ability to preserve constant states for any mesh motion. This is related to the conservation of cell volumes in relation to the mesh motion. In our scheme, if we start with a constant state , then the predictor is also constant in the spacetime interval, i.e., . The spacetime terms in (3.2) are polynomials with degree in space and degree one in time and these are exactly integrated by the chosen quadrature rule. The flux terms at cell boundaries in (3.2) are of degree one in time and these are also exactly integrated. Hence the scheme (3.2) can be written as
where and for . Due to the constant predictor and by consistency of the numerical flux
Moreover, for
where we have used the property that is an affine function of and are orthogonal. This implies that for . For , we get
and since , we obtain which implies that .
8 Grid coarsening and refinement
The size of the cells can change considerably during the time evolution process due to the near Lagrangian movement of the cell boundaries. Near shocks, the cells will be compressed to smaller sizes which will reduce the allowable time step since a CFL condition has to be satisfied. In some regions, e.g., inside expansion fans, the cell size can increase considerably which may lead to loss of accuracy. In order to avoid too small or too large cells from occurring in the grid, we implement cell merging and refinement into our scheme. If a cell becomes smaller than some specified size , then it is merged with one of its neighbouring cells and the solution is transferred from the two cells to the new cell by performing an projection. If a cell size becomes larger than some specified size , then this cell is refined into two cells by division and the solution is again transferred by projection. The use of projection for solution transfer ensures the conservation of mass, momentum and energy and preserves the accuracy in smooth regions. We also ensure that the cell sizes do not change drastically between neighbouring cells. To keep a track of refinement of cells, each cell is assigned an initial level equal to . The daughter cells created during refinement are assigned a level incremented from the parent cell, while the coarsened cells are assigned a level decremented from the parent cell.
The algorithm for refinement and coarsening is carried out in three sweeps over all the active cells. In the first sweep, we mark the cells for refinement or coarsening based on their size and the level of neighboring cells. Cells are marked for coarsening if the size is less than a prespecified minimum size. They are marked for refinement if either the size of the cell is larger than the maximum size or if the level of the cell is less than the level of the neighboring cells. If none of the conditions are satisfied, the cells are marked for no change. In the second sweep, a cell is marked for refinement if both the neighboring cells are marked for refinement. A cell is also marked for refinement if the size of the cell is larger than twice the size of either of the neighboring cells, and is also larger than twice the minimum size. The last condition is inserted in order to prevent a cell being alternately marked for refinement and coarsening in consecutive adaptation cycles. In the third and final sweep, we again mark cells for refinement if both the neighboring cells are marked for refinement. Further, we ensure that a cell marked for refinement does not have a neighboring cell marked for a coarsening, since this can lead to an inconsistent mesh.
9 Numerical results
The numerical tests are performed with polynomials of degree one, two and three, together with the linear Taylor expansion, two stage CERK and four stage CERK, respectively, for the computation of the predictor. For the quadrature in time, we use the midpoint rule, two and three point GaussLegendre quadrature, respectively. The time step is chosen using the CFL condition,
where is given by equation (10), and the factor comes from linear stability analysis cockburn1989 . In most of the computations we use unless stated otherwise. We observe that the results using average or linearized Riemann velocity are quite similar. We use the average velocity for most of the results and show the comparison between the two velocities for some results. The main steps in the algorithm within one time step are as follows.

Choose mesh velocity

Choose time step

Compute the predictor

Update solution to the next time level

Apply TVD/TVB limiter on

Apply positivity limiter on from Zhang:2010:PHO:1864819.1865066 .

Perform grid refinement/coarsening
In all the solution plots given below, symbols denote the cell average value.
9.1 Order of accuracy
We study the convergence rate of the schemes by applying them to a problem with a known smooth solution. The initial condition is taken as
whose exact solution is , , . The initial domain is and the final time is units. The results are presented using Rusanov and HLLC numerical fluxes. The norm of the error in density are shown in table (1), (3) for the static mesh and in table (2), (4) for the moving mesh. In each case, we see that the error behaves as which is the optimal rate we can expect for smooth solutions. In table (5), we show that the ALE DG methods preserves its higher order in presence of a limiter.
The mesh velocity is constant since the fluid velocity is constant. In order to study the effect of perturbations in mesh velocity, we add a random perturbation to each mesh velocity, where
is a uniform random variable in
and and a sample velocity distribution is shown in figure (5). Note that this randomization is performed in each time step with different random variables drawn for each face. For moving mesh, there is no unique cell size, and the convergence rate is computed based on initial mesh spacing which is inversely proportional to the number of cells. From table (6) which shows results using HLLC flux, we again observe that the error reduces at the optional rate of even when the mesh velocity is not very smooth.Error  Rate  Error  Rate  Error  Rate  

100  4.370E02    3.498E03    3.883E04   
200  6.611E03  2.725  4.766E04  2.876  1.620E05  4.583 
400  1.332E03  2.518  6.415E05  2.885  9.376E07  4.347 
800  3.151E04  2.372  8.246E06  2.910  5.763E08  4.239 
1600  7.846E05  2.280  1.031E06  2.932  3.595E09  4.180 
Error  Rate  Error  Rate  Error  Rate  

100  2.331E02    3.979E03    8.633E04   
200  6.139E03  1.9250  4.058E04  3.294  1.185E05  6.186 
400  1.406E03  2.0258  5.250E05  3.122  7.079E07  5.126 
800  3.375E04  2.0366  6.626E06  3.077  4.340E08  4.760 
1600  8.278E05  2.0344  8.304E07  3.057  2.689E09  4.573 
Error  Rate  Error  Rate  Error  Rate  

100  4.582E02  3.952E03  3.464E04  
200  9.611E03  2.253  4.048E04  3.287  2.058E05  4.073 
400  2.052E03  2.240  4.640E05  3.206  1.287E06  4.036 
800  4.803E04  2.192  5.623E06  3.152  8.061E08  4.023 
1600  1.184E04  2.149  6.929E07  3.119  5.050E09  4.016 
Error  Rate  Error  Rate  Error  Rate  

100  1.590E02  1.626E03  1.962E04  
200  4.042E03  1.977  2.072E04  2.972  1.269E05  3.950 
400  1.014E03  1.985  2.605E05  2.982  7.983E07  3.971 
800  2.538E04  1.990  3.261E06  2.988  4.997E08  3.980 
1600  6.349E05  1.992  4.077E07  2.991  3.124E09  3.985 
Error  Rate  Error  Rate  

100  2.053E02    2.277E03   
200  4.312E03  2.251  3.425E04  2.732 
400  1.031E03  2.064  4.565E05  2.907 
800  2.550E04  2.015  5.812E06  2.973 
1600  6.356E05  2.004  7.315E07  2.990 
Error  Rate  Error  Rate  Error  Rate  

100  1.735E02    1.798E03    2.351E04   
200  4.179E03  2.051  2.848E04  2.676  1.416E05  4.069 
400  1.054E03  2.035  4.301E05  2.703  8.578E07  4.041 
800  2.615E04  1.943  6.012E06  2.838  5.476E08  3.958 
1600  7.279E05  1.852  8.000E07  2.909  3.505E09  3.966 
9.2 Smooth Test Case with NonConstant Velocity
We also test the accuracy of our schemes on a isentropic problem with smooth solutions. The test case The initial conditions are given by
Comments
There are no comments yet.