This paper considers the numerical solutions of hyperbolic systems of conservation laws on discontinuous grids. The discontinuity in the grid is often due to the overlapping of different mesh systems, or to the adaptive mesh refinement. In either case, if the discontinuities in solutions are genuine nonlinear, some oscillations will be conserved around the grid interface for nonlinear systems.
Computations on discontinue grids are becoming more common for two reasons. First, many multidimensional problems of practice interest involve complex geometry, and in general it is not sufficient to be able solve hyperbolic equations on a uniform Cartesian grid in a rectangular domain. As the configurations that can be modeled become more complex, so does the grid generation problem. Generally, it is very difficult to generate one smooth body-fitted grid to cover the whole complex domain. To simplify this procedure, it is accustomed to use multigrid to fit the complex domain. Each part of the domain will still have a smooth grid, but now the component grids will in general overlap rather than in an irregular fashion.
The second reason of these discontinue grids comes from the use of adaptive methods     , where a grid is abruptly refined in order to gain high- resolution solution around areas of large solution variation. In this approach, many fine grids are distributed on those regions where high resolution is required. Efficient implementation of the adaptive strategy can increase the accuracy of the numerical approximations and meanwhile decrease the computational cost.
The main objective of this paper is to investigate postshock oscillations around the grid interface. Some linear and nonlinear equations are tested with LLF scheme, and this kind of oscillations only appears in nonlinear equation system. In addition, Godunov scheme behaves similarly to the LLF scheme, and so we don’t give the numerical results of Godunov scheme for saving space. We incline to believe that the postshock oscillations are an inevitable feature of shocks captured by currently frequently-used methods. The postshock oscillations due to discontinue grids are different from common nonphysical oscillations caused by high-order interpolation, but it is very similar to that appear behind the slow shocks   . The amplitude and concave-convex shape of oscillations are determined by the ratio between coarse- and fine-grid length. From the point of Riemann invariant, these oscillations have similar feature with rarefaction. That is Riemann invariant has the same value at all points on the corresponding oscillation.
The rest of the paper is organized as follows. In Section 2, we briefly review some model problems and give the mesh distribution. In Section 3, the LLF scheme is described. The numerical results and the observed behavior are demonstrated in Section 4. However, we have yet found the reason which leads to the postshock oscillations observed in the paper.
2. Equations and the Grid Distribution
The postshock oscillations are investigated under the following five sets of equations: the advection equation,
the inviscid Burgers equation,
the shallow water equation,
where is the depth of the fluid, is the horizontal velocity and is acceleration of gravity; and full Euler equations for an ideal gas with constant specific heats,
where are the density, velocity, total energy and pressure, respectively. The system of Euler equations is closed by the equation of state for an ideal polytropic gas
here is the ratio of specific heat.
The last example is 2D Euler equations,
where the above system is closed by the equation of state,
For ease of presentation, we consider the case of mesh refinement by a factor . Figure 1 shows a grid consists of two uniform subgrids and the factor . Suppose the problem is on the physical domain , which is subdivided into cells with and , so that
and . Apparently, the interface between the coarse and fine grids is located at .
3. Numerical Methods
The equations presented in the previous section will be solved using first order accurate LLF  method. The scheme for solving Equations (1)-(5) can be written in the form
where is the numerical flux function at the interface between cells and at time level . Depending on the choice of the numerical flux formula, Equation (6) is referred to as LLF scheme.
LLF scheme is the improved version of the classical Lax-Friedrichs scheme by replacing the value by a locally determined value,
Figure 1. Interface between the coarse and fine grids.
4. Numerical Result
4.1. Advection Equation
Firstly, advection Equation (1) is considered on a non-uniform grid which is composed of two uniform subgrids with initial condition
The ratio between the coarse and fine grids is , and the interface in the physics domain is located at . The boundary of domain , , and the number of cell , . Obviously, the discontinuity will arrive the interface when and wholly pass through it when .
For the linear advection equation, LLF scheme reduces to upwind scheme. The solutions to advection equation can be seen in Figure 2. From this figure, we find that the solution has monotone profile as expected and no oscillation can be observed around .
4.2. Inviscid Burgers Equation
The second example is the inviscid Burgers Equation (2) subject to the initial data
The solution propagates to the right with shock speed . Figure 3 shows the solution at , and the shocks pass through the interface at . The zoomed version around is shown in Figure 3. As in the previous problem, LLF scheme shows monotone solutions which are similar to the solutions on the uniform grids.
4.3. Shallow Water Equations
Consider the shallow water equations (3) with the piecewise-constant initial conditions
This is a special case of the Riemann problem in which , and is called the dam-break problem because it models what happens if a dam separating two levels of water bursts at time . This is the shallow water equivalent of the shock tube problem of gas dynamics.
Note that the Jacobian matrix of the shallow water equations is
The eigenvalues of are
Figure 2. Left: numerical solution at with LLF scheme. Right: zoomed region around the interface between the coarse and fine grids.
Figure 3. Left: numerical solution of Burgers equation at with LLF scheme. Right: zoomed region around the interface between the coarse and fine grids.
with the corresponding eigenvectors
Firstly, a uniform grid is considered. Figure 4 shows the structure of solution of this Riemann problem. Note that the shallow water equations are a system of two equations, and so the Riemann solutions generally contain two waves. For this initial-value problem (3) (8), these always consist of one 1-rarefaction wave on the left and one 2-shock on the right with shock speed . In this case, we find that the solutions behave well and no postshock oscillation is introduced.
For the shallow water equations, Leveque  gives detailed formula of Hugoniot
Figure 4. Left: solution of the dam-break Riemann problem for the shallow water equations with on uniform grids; Right: zoomed version around .
loci and integral curves. The Hugoniot curves correspond to the 1-shock and 2-shock are
The integral curves correspond to 1-rarefaction and 2-rarefaction waves are
As shown in previous paragraph, the solution of this initial-value problem contains both shock and rarefaction. So, Hugoniot and integral curves originate from the left state and right state , and the intersection of them is intermediate state . Figure 5 shows that numerical solution by LLF scheme almost lies on the Hugoniot and integral curves.
In Figure 6, physics domain is divided into two subdomain and , and each subdomain is distributed some sort of uniform grids. Denote the length of the cells on interval and are and respectively. Figure 6 shows the solution on nonuniform grids with at . Amplifying the figure around interface at , an obvious oscillation can be observed behind the shock (on the left of the interface). Also, we find that the shock is in the correct location and the shock speed remain unchanged
Figure 5. Left: Hugoniot loci and integral curve of the eigenvector for the shallow water equations in the state space ; Right: zoomed version.
Figure 6. Solution of the dam-break Riemann problem for the shallow water equations with nonuniform grids which consist of two uniform grids with .
after it pass through the interface. The oscillation is the wave that arose from the initial discontinuity at when the shocks pass the interface.
A interesting phenomenon can be observed in Figure 7. The hump (oscillation) in Figure 6 changes into dip if we shift from 2 into 0.5. Note that the hump and dip both belong to acoustic wave moving at speed , where .
Note that, for the shallow water equations, the function of 1-Riemann invariant and 2-Riemann invariant are
Figure 8 shows the Riemann invariants of the equations. Generally, the function of 1-Riemann invariant takes the same value at all points on the integral curve of . Similarly, the value of 2-Riemann invariant is constant along any integral curve of . We find that the acoustic wave (oscillation behind the shock) takes the same value as 1-rarefaction wave for 1-Riemann invariant.
Figure 9 shows Hugoniot, integral curves and numerical orbit in state space on nonuniform grids. Note that the numerical orbit deviating the integral curve in Figure 9 is due to the departure between numerical solution and exact solution at the left of the shock. In addition, a redundant tip appears in this figure which is different from the distribution of numerical orbit in Figure 5.
4.4. 1D Euler Equations
Euler equations of gas dynamics is discussed in this section. We test the postshock oscillations with initial data
This is a well known test problem proposed by Sod  . The physical domain
Figure 7. Solution of the dam-break Riemann problem for the shallow water equations with nonuniform grids which consist of two uniform grids with .
Figure 8. The numerical orbit of the Riemann invariants for the shallow water equations. Left: 1-Riemann invariant; Right: 2-Riemann invariant.
Figure 9. Left: Hugoniot loci and integral curve of the eigenvector for the shallow water equations in the state space ; Right: zoomed version.
Figure 10. Numerical solution of Sod problem on nonuniform grids at .
is also divided into two uniform subdomain and .
Except for the first and third characteristic fields are genuinely nonlinear which have similar behavior to the two characteristic fields in the shallow water equations, the second field is linearly degeneration. So, the solution to Riemann problem for Euler problem typically has two nonlinear waves and a contact discontinuity.
Figure 10 shows the numerical solution on nonuniform grids with . The shock has already passed through the interface at , and two humps arise near . This is different from what we observed in shallow water equations. In that system, there is just one wave arises when the shock passes through the interface between the coarse and fine grids, but here two waves arise for Euler equations. As mentioned previously, the Riemann problem for shallow water and Euler equations generally consist of two and three waves respectively. So, the phenomena that postshock oscillations arose when the shock through the interface is similar to the solution of Riemann problem.
From Figure 10, we find that two waves propagate with different direction. Note that these two waves are different from some nonphysical oscillations. The left wave propagate to the left with the speed , and the left wave propagate to the other direction with the speed , where asterisk denotes intermediate state and . In addition, the amplitude of vibration is decided by the ratio . The bigger ratio between the coarse and fine grids, the bigger amplitude the postshock will have.
For this initial value problem, all of the Riemann invariants for a polytropic ideal gas are summarized below:
Figure 11 shows the 1-Riemann invariants for this initial value problem. We
Figure 11. 1-Riemann invariants for Euler equations with initial data (9).
Figure 12. 2-Riemann invariants.
find that the 1-Riemann invariants are constant across the rarefaction and the left wave which generate from , and only one oscillation was left. Similarly, Figure 12 shows the 2-Riemann invariants which are constant across the second oscillation at the right of . And 3-Riemann invariant is given in Figure 13. From these three figures, the oscillations behind the shock behave similarly to the rarefaction for Riemann invariants.
4.5. 2D Euler Equations
The last example is 2D Euler equations (5), with initial data
which corresponds to the case of left forward shock, right backward shock, upper backward shock, and lower forward shock. We refer the readers to    for details. This example can be used to test the postshock oscillations as shock across the interface.
Figure 14 shows the cross-sectional solution at . One solution on uniform grids and another solution on nonuniform grids are shown in this figure. The obvious difference between these two solutions is that two additional humps appear in the solution under nonuniform grids. This phenomenon is similar to what appear in the previous nonlinear system.
The postshock oscillations appeared in the nonlinear system has been elucidated numerically. We found that these oscillations would not appear in the nonlinear scalar equation and linear equations system. Although all the numerical results were computed with LLF method, these oscillations also appeared with other schemes, such as Godunov, Roe schemes and so on. So the postshock oscillations
Figure 13. 3-Riemann invariant.
Figure 14. Cross section at . Left: the amplified figure.
have nothing to do with the specific methods which are used to compute solution on nonuniform grids. In addition, this kind of oscillations is different from those nonphysical oscillations arose near the shock, and they are similar to those emerge behind the slow-moving shocks which can be found in   .
This work was supported by the National Natural Science Foundation of China (No. 11501238), Natural Science Foundation of Guangdong Province (No. 2012A030313119) and Supported by the Major Project Foundation of Guangdong Province Education Department (No. 2014KZDXM070).
 Berger, M.J. and Oliger, J. (1984) Adaptive Mesh Refinement for Hyperbolic Partial Differential Equations. Journal of Computational Physics, 53, 484-512.
 Huang, W., Ren, Y. and Russell, R.D. (1994) Moving Mesh Partial Differential Equations (MMPDES) Based on the Equidistribution Principle. SIAM Journal on Numerical Analysis, 31, 709-730.
 Stockie, J.M., Mackenzie, J.A. and Russell, R.D. (2001) A Moving Mesh Method for One-Dimensional Hyperbolic Conservation Laws. SIAM Journal on Scientific Computing, 22, 1791-1813.
 Sod, G.A. (1978) A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws. Journal of Computational Physics, 107, 1-31.
 Lax, P.D. and Liu, X.D. (1998) Solutions of Two-Dimensional Riemann Problems of Gas Dynamics by Positive Schemes. SIAM Journal on Scientific Computing, 19, 319-340.
 Schulz-Rinne, C.W., Collins, J.P. and Glaz, H.M. (1993) Numerical Solution of the Riemann Problems for Two-Dimensional Gas Dynamics. SIAM Journal on Scientific Computing, 14, 1394-1414.
 Tang, H.Z. and Tang, T. (2003) Adaptive Mesh Methods for One-and Two-Dimensional Hyperbolic Conservation Laws. SIAM Journal on Numerical Analysis, 41, 487-515.