Conservation error for 1-D Euler with WENO

kiliancfd · July 22, 2014, 19:01

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.

praveen · July 23, 2014, 06:46

We can write the Roe flux as

$F_{j+1/2} = \frac{1}{2}[F(u_{j+1/2}^L) + F(u_{j+1/2}^R)] - \frac{1}{2}D(u_{j+1/2}^L, u_{j+1/2}^R) (u_{j+1/2}^R - u_{j+1/2}^L)$

where the matrix D involves the Roe averages. You are saying this does not lead to full convergence, but if you use the following form

$F_{j+1/2} = \frac{1}{2}[F(u_{j+1/2}^L) + F(u_{j+1/2}^R)] - \frac{1}{2}D(u_{j}, u_{j+1}) (u_{j+1/2}^R - u_{j+1/2}^L)$

you get full convergence and exact solution. Does this describe your situation ?

Note that the second form is also consistent and many people use it.

kiliancfd · July 23, 2014, 13:54

I would not say that using the second form gives full convergence, since the solution never changes at any time step, but otherwise that does describe my situation.

With the first formulation the momentum in the solutions is not conserved. With the second formulation the computations "converge" to the exact solution but only because that is the initial condition I provided. It should converge to the solution of the discretized problem which has a smeared shock (since the Lax-Friedrichs splitting introduces some dissipation) and constant momentum in the whole interval.

praveen · July 23, 2014, 23:51

Since you use a finite volume method total momentum is conserved. Why do you say it is not conserved ?

In some cases, you can get exact solution with riemann solvers. Maybe you are in that situation.

kiliancfd · July 24, 2014, 00:49

The scheme does conserve momentum as you say, but the result does not show the momentum being constant across the shock. That is what I mean by saying that the momentum is not conserved.

I've attached a picture of my results in which you can see that the momentum (black line) peaks around the shock (marked by the blue line which is the density). My understanding is that this peak should not exist, but if I use the first formulation of the Roe flux then it always appears.

praveen · July 24, 2014, 03:27

You may not get the momentum to be constant but the fluxes should be nearly constant. Can you check if the spike decreases with grid refinement. If it does then you dont need to worry about it.

kiliancfd · July 24, 2014, 17:26

Unfortunately the peak does not become smaller on finer grids. However, it does shrink when I decrease the upstream Mach number.

In my code I do a characteristic decomposition that treats the eigenvector matrix as locally constant. Could this spike arise because the eigenvectors are not constant across the shock? The impression I have is that this local constancy is always assumed.

July 22, 2014, 19:01	Conservation error for 1-D Euler with WENO	#1
kiliancfd New Member Kilian Join Date: Jul 2014 Location: United States Posts: 4 Rep Power: 11	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.

July 23, 2014, 06:46		#2
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 342 Blog Entries: 6 Rep Power: 18	We can write the Roe flux as $F_{j+1/2} = \frac{1}{2}[F(u_{j+1/2}^L) + F(u_{j+1/2}^R)] - \frac{1}{2}D(u_{j+1/2}^L, u_{j+1/2}^R) (u_{j+1/2}^R - u_{j+1/2}^L)$ where the matrix D involves the Roe averages. You are saying this does not lead to full convergence, but if you use the following form $F_{j+1/2} = \frac{1}{2}[F(u_{j+1/2}^L) + F(u_{j+1/2}^R)] - \frac{1}{2}D(u_{j}, u_{j+1}) (u_{j+1/2}^R - u_{j+1/2}^L)$ you get full convergence and exact solution. Does this describe your situation ? Note that the second form is also consistent and many people use it. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

July 23, 2014, 23:51		#4
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 342 Blog Entries: 6 Rep Power: 18	Since you use a finite volume method total momentum is conserved. Why do you say it is not conserved ? In some cases, you can get exact solution with riemann solvers. Maybe you are in that situation. FMDenaro likes this. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

July 24, 2014, 03:27		#6
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 342 Blog Entries: 6 Rep Power: 18	You may not get the momentum to be constant but the fluxes should be nearly constant. Can you check if the spike decreases with grid refinement. If it does then you dont need to worry about it. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
WENO code for Euler equations	buaalzr	Main CFD Forum	19	October 16, 2017 16:32
Euler Equations, Sod shock tube & conservation	Antigravity	System Analysis	1	July 17, 2014 13:30
Question about WENO limiter	DJChoi	Main CFD Forum	0	August 21, 2012 10:03
euler euler multiphase modelling	mdsvjpk	FLUENT	0	June 14, 2011 11:50
convergence of Euler Equations with WENO	Shyam	Main CFD Forum	0	June 1, 2007 00:47

July 23, 2014, 13:54		#3
kiliancfd New Member Kilian Join Date: Jul 2014 Location: United States Posts: 4 Rep Power: 11	I would not say that using the second form gives full convergence, since the solution never changes at any time step, but otherwise that does describe my situation. With the first formulation the momentum in the solutions is not conserved. With the second formulation the computations "converge" to the exact solution but only because that is the initial condition I provided. It should converge to the solution of the discretized problem which has a smeared shock (since the Lax-Friedrichs splitting introduces some dissipation) and constant momentum in the whole interval.

July 24, 2014, 17:26		#7
kiliancfd New Member Kilian Join Date: Jul 2014 Location: United States Posts: 4 Rep Power: 11	Unfortunately the peak does not become smaller on finer grids. However, it does shrink when I decrease the upstream Mach number. In my code I do a characteristic decomposition that treats the eigenvector matrix as locally constant. Could this spike arise because the eigenvectors are not constant across the shock? The impression I have is that this local constancy is always assumed.