Woodward-Collela blast wave problem with Finite Volume WENO

We consider the one-dimensional Euler equations for inviscid, compressible flow: \[\begin{align*} \rho_t + (\rho u)_x &= 0 \\ (\rho u)_t + (\rho u^2 + p)_x &= 0 \\ E_t + (u (E + p) )_x &= 0, \end{align*}\] where \(\rho\) is the density, \(u\) the velocity, \(E\) the total energy and \(p\) the pressure. The fluid is an ideal gas, thus the pressure given by \(p=\rho(\gamma - 1)e\), where \(e\) is internal energy and \(\gamma = 1.4\). The spatial domain is the unit interval and for initial conditions we choose \[\begin{align*} \rho(x,0) = 1, \quad u(x,0) = 0, \quad \text{ and } p(x,0) = \begin{cases} 1000 & \mbox{if } 0 \le x < 0.1 \\ 0.01 & \mbox{if } 0.1 \le x < 0.9 \\ 100 & \mbox{if } 0.9 \le x \le 1. \end{cases} \end{align*}\] The initial conditions consist of two parallel, planar flow discontinuities. Reflecting boundary conditions are used and the two discontinuities yield strong shock waves and contact discontinuities going inwards the domain, whereas the rarefactions going outwards are reflected and interact with the other waves.

In CLAWPACK, a script that runs the Woodward-Collela blast problem can be found at clawpack/pyclaw/examples/euler_1d. Stable numerical solutions for this problem are found with second order TVD discretization coupled with strong stability preserving Runge Kutta methods (SSPRK). We would like to use finite volume WENO discretizations instead. Several attempts have been made, but all fail to find a solution after the two shocks collide. Note that in the literature finite difference WENO with SSPRK methods have been used, but as far we are concerned no finite volume WENO has been used.

All attempts listed below and can be found at https://github.com/hadjimy/pyclaw/tree/char_decomp_blast_prb and can be reproduced by running woodward_collela_blast.py in the char_decomp_weno repository.

1. No characteristic decomposition:

This routine implements Chi-Wang Shu’s code for fifth order finite volume WENO reconstruction over the conservative variables. The reconstruction is performed component-wise and it is the default option in CLAWPACK when no characteristic decomposition is used. (see subroutine weno5 of reconstruct.f90; lines 158–228)

2. No characteristic decomposition - reconstruction over pressure:

This routine implements fifth order finite volume WENO reconstruction on the velocity, momentum and pressure instead on the conservative variables. Again no characteristic decomposition is used and after reconstruction of pressure the energy is computed so that reconstruction of conservative variables is completed (in a implicit way). (see subroutine weno5_pressure of reconstruct.f90; lines 232–354)

3. Reconstruction over primitive variables:

Change from conservative variables to primitive variables and perform fifth order finite volume WENO reconstruction over primitive variables . No characteristic decomposition is used and after reconstruction we change back to conservative variables. (see subroutine weno5_primitive of reconstruct.f90; lines 358–451)

4. Characteristic decomposition over primitive variables:

The routine changes from conservative variables to primitive variables. Left and right eigenvector matrices are computed over the primitive variables and fifth order finite volume WENO characteristic reconstruction is performed. Then it projects to the physical space and changes back to conservative variables. (see subroutine weno5_char_primitive of reconstruct.f90; lines 455–571)

5. Characteristic decomposition over cell averages:

This routine does projection on the characteristic space, then performs WENO reconstruction andprojects back to physical space. Left and right eigenvector matrices at each interface are computed over the Roe averages. Note that characteristic projections are computed over cell averages instead of differences of the cell averages. This is not the correct way for characteristic-wise WENO and adds difussion to the problem (the correct is in weno5\_char). (see subroutine weno5_char_cell_avg of reconstruct.f90; lines 575–675)

6. Characteristic decomposition over over differences of the cell averages:

This routine implements projection on the characteristic space and performs WENO reconstruction based on Chi-Wang Shu’s code. Then it projects back to physical space. Left and right eigenvector matrices at each interface are computed over the Roe averages. Note that characteristic projections are computed over difference of cell averages. This is the default routine for characteristic decomposition in CLAWPACK. (see subroutine weno5_char of reconstruct.f90; lines 679–775)

7. Characteristic decomposition over over differences of the cell averages (expanded reconstruction):

This routine implements projection on the characteristic space and performs WENO reconstruction in an equivalent way to Chi-Wang Shu’s code. The difference from weno5\_char is on the linear part of the reconstruction which is stenil-independent and can be computed separately from the characteristic reconstruction.
Then the routine projects back to physical space. Left and right eigenvector matrices at each interface are computed over the Roe averages and characteristic projections are computed over difference of cell averages. (see subroutine weno5_char_clean of reconstruct.f90; lines 779–904)

Remark 1: When no characteristic decomposition is used, solution exhibits small oscillation (before the two shocks collide), hence it is important to note that characteristic decomposition is essential for an accurate solution.

Remark 2: When characteristic decomposition is used the left and right eigenvector matrices at each interface can be computed over cell averages or over the arithmetic mean of cell averages or over the Roe mean of cell averages. The correct option is to use Roe averages of the conservative (or primitive) variables since a Roe solver is used. However, the other options seems to have no visible effects.