[kiliancfd]

I am using a WENO method to solve the 1-D Euler equations with the ideal gas equation of state. As a test problem I am using a stationary shock where the initial conditions are those for a normal shock with an upstream Mach number of 2. The shock smears over a few grid points but otherwise retains the piecewise constant values from the initial condition, which is as expected, but the residuals (i.e., the average absolute value of the rate of change of density) do not decrease to machine zero. Instead, the residual remains at around 10^-6 indefinitely.

I have found that the density changes only at the four points that resolve the shock. At those same points, the momentum deviates from the constant value that it should have across the entire domain and I believe this failure of conservation is responsible for the residuals' hanging. I do not, however, see why that violation arises in the first place. It appears after the very first time step and keeps the same shape throughout the computations.

Has someone here encountered a similar issue before, or have some idea what the cause might be?


A couple of technical details to explain my situation more fully:

1) The basic outline of the reconstruction procedure is as follows:
1) For each face use the cell-average values to either side to calculate a Roe-averaged state, from which the eigenvalues and eigenvectors of the flux Jacobian at the interface are computed and used in a characteristic decomposition.
2) Use a WENO scheme to obtain left- and right-biased reconstructions of the characteristic variables at the current face.
3) From these biased estimates, compute a new Roe-averaged state vector.
4) Compute the eigenvalues of the corresponding flux Jacobian and use them to upwind appropriately. I use the Roe-fixed approach, which incorporates the local Lax-Friedrichs flux splitting as an entropy fix for the Roe scheme (specifically, Eq. 5.7 from DOI: 10.1137/110857659).

Remark: Intuitively, it makes sense to me that the second Roe-average should use the two reconstructions rather than the known cell averages. However, when I modify that step to use the cell averages the problem I described does not occur. In fact nothing at all occurs; the solution remains precisely equal to the initial condition at all times. So neither approach really works well.

2) I am using the smoothness indicators and weights proposed by Jiang & Shu, with epsilon = 10^-7, but the same trouble happens for other values of epsilon and other weight/indicator combinations. I have tried those of Borges et al. and Henrick, Aslam & Powers with the same result.

Thoughts or suggestions? Thank you for your time.