Implementation of the Roe and Riemann Exact Solver

Godunov-type schemes and all its higher-order generalisations require the solution of a Riemann problem for every cell boundary at each time-step. The exact iterative Riemann solver of Gottlieb and Groth [13] typically needs the computation of an initial guess solution for the Newton iterative method, the identification of the pattern of the solution and finally some Newton iterations. It is clearly a computationally expensive procedure. Moreover, most of the structure of the resulting solution is lost in the Godunov's method, due to the final cell-average operation on each grid-cell. This suggests that good numerical results could be obtained by calculating an appropriate approximate solution instead of the exact one.

Roe's algorithm [14] solves exactly a linearized problem, instead of looking for an iterative solution of the exact original Riemann problem. The approximate solver proposed by Roe is much less expensive in terms of computational effort than the exact one, because the exact solution of a linear Riemann problem can be more easily built. This is the reason why we do not need Burgers equation to figure out the exact solution of the triangular profile (it is not straightforward tough).

When the energy of a flow is highly kinetic, a Godunov-type conservative numerical scheme based on a veneration of the Riemann problem may fail by predicting non-physical states with negative density or internal energy. This failure may happen because in flows where the dominant energy mode is kinetic the resulting internal energy may be negative. When the jump in the shock produces a null zero shock speed, i.e. vanishing, we expect differences. In these cases a Roe linearization will fail and produce different fluxes than the Riemann exact solver.

In addition to characteristics crossing and a shock forming, another possible way the method of characteristics can break down, and a discontinuity can form, is if the characteristics on both sides of the discontinuity emanate from it, rather than go into it. Then the discontinuity is called a rarefaction wave.

For reasons of numerical efficiency and, particularly in multi-dimensions, the exact solution of the Riemann problem is frequently avoided and linearized (approximate) Riemann solvers are preferred. These solvers are based on the exact solution of Riemann problems corresponding to a linearized version of the original system of equations. After extensive experimentation, the results achieved with approximate Riemann solvers are comparable to those obtained with the exact solver.

We modified calc_flux.f so the flux is calculated now via:

     f(1) = 0.5d0*q**2

And to implement the exact Riemann solver we wrote a subroutine named calc_Riemann_exact.f which calculates the flux using the following algorithm (Ian Hawke's thesis page 14 [12]),

$${\bf f(q)} = \left\{ \begin{array}{rl} \textbf{if $q_L < q... ...c{1}{2}(q_L+q_R) < 0$} \end{array} \right. \end{array}\right.$$

This new subroutine is called by get_flux.f.