Exact Riemann Solver (no vacuum)

The Riemann problem is the solution of a system of hyperbolic conservation laws (in our case the Euler equations). This solver assumes two one-dimensional ideal gases separated at the origin x=0, with initial densities, velocities, and pressures given as input. It gives a solution to the set of equations

\rho + \rho u = 0$ $\rho u + \rho u^2 + p = 0$ $E + u(E + p) = 0$

with initial conditions

\rho(x,0) = \begin{cases} \rho_L, & x < 0 \\ \rho_R, & x > 0 \end{cases}$ $\rho(x,0)u(x,0) = \begin{cases} \rho_L u_L, & x < 0 \\ \rho_R u_R, & x > 0 \end{cases}$ $E(x,0) = \begin{cases} E_L, & x < 0 \\ E_R, & x > 0 \end{cases}$

The solutions to this IVP at later times are split into four regions - to the left, the initial left state, unperturbed by interaction with the right; between the left and center-left regions, there will either be a rarefaction fan or a shock (assumed ideal), across which the new equilibrium left state holds; the left and right fluids will be separated by a contact discontinuity, across which the new equilibrium right state holds (note that the right and left equilibrium states differ only in their densities); the center-right region will then be separated by a rarefaction fan or (again ideal) shock from the rightmost region, which is still in the initial right state.

The program reads input from a file called problem.data, and stores many of the relevant parameters in global variables.

NAMELIST/ProblemData/ rhoL,uL,pL,rhoR,uR,pR,gamma,tol,starttime,finaltime,frames,xL,xR,res

inquire(file='problem.data', exist=ex)
if (ex) then
 OPEN(UNIT=30, FILE='problem.data', STATUS="OLD")
 READ(30,NML=ProblemData)                               ! read into variables in namelist
 CLOSE(30)
else
 print *, 'Missing problem.data. Exiting.'
 stop           ! exit if missing
end if

The state vectors W_L and W_R are created with overloaded constructors that also calculate useful parameters A and B and the speed of sound c for use later in the program.

These state vectors are then passed to the Newton-Raphson solver, which calls a subroutine to guess a good value for p₀ (currently not very effective - just uses tol) and then iterates until the tolerance for change in pressure between steps is met.

f_reg(p_*) = (p_* - p_reg) [\frac{2}{(\gamma + 1) \rho_reg (p_* + \frac{\gamma - 1}{\gamma + 1} p_reg)}]$

See below for derivation of function f.

From the shock jump conditions

\rho_L \hat{u_L} = \rho_{*L} \hat{u_*}$ $\rho_L \hat{u_L}^2 + p_L = \rho_{*L} \hat{u_*}^2 + p_*$ $\hat{u_L}(\hat{E_L} + p_L) = \hat{u_*}(\hat{E_{*L}} + p_*)$

and the energy and EOS of an ideal gas

e = \frac{p}{(\gamma - 1) \rho}$ $p = \rho R T$

where \hat{u} = u_L - S_L is the velocity of the fluid in the frame of the shock, the function f_L(p_*) = u_L - u_*,

f_L(p_*) = (p_* - p_L) [\frac{2}{(\gamma + 1) \rho_L (p_* + \frac{\gamma - 1}{\gamma + 1} p_L)}]$

can be derived to connect the left state with the pressure in the central, interaction region in the case of a shock. Similarly, in the case of a rarefaction, instead of the shock jump conditions, the isentropic equation

\frac{p_*}{\rho_*} = \frac{p_L}{\rho_L}$

and the Generalized Riemann Invariant

u_L + \frac{2 c_L}{\gamma - 1} = u_* + \frac{2 c_{*L}}{\gamma - 1}$

where c = \sqrt{\frac{\gamma p}{\rho}} is the speed of sound in the region, can be used to construct

f_L(p_*) = \frac{2 c_L}{\gamma - 1}[(\frac{p_*}{p_L})^{\frac{\gamma - 1}{2 \gamma}} - 1]$

The same analysis can be performed on the right side of the system to yield the same results.

The Riemann problem can then be solved by finding the value p = p_* such that

f_L(p_*) + f_R(p_*) + u_R - u_L = 0$

The function and its first and second derivatives are well-behaved and monotonic, so we can use a simple Newton-Raphson iterative method, derived from a Taylor series expansion of f(p_*):

If our initial guess is p_k and f(p_k) \ne 0, we want to change p_k so that f(p_k + \delta) = f(p_{k+1}) = 0; the Taylor expansion is given by

f(p_k + \delta) = f(p_k) + \delta f'(p_k) + O(\delta^2)$

Dropping second order and higher terms, since we have f(p_k + \delta) = 0, we know f(p_k) + \delta f'(p_k) = 0, and solving for \delta,

\delta = -\frac{f(p_k)}{f'(p_k)}$

and

p_{k + 1} = p_k + \delta = p_k - \frac{f(p_k)}{f'(p_k)}$

The iteration procedure is terminated at the term such that

|\frac{2 (p_{k + 1} - p_k)}{p_{k + 1} + p_k}| < Tol$

where Tol is the (generally small) tolerance passed to the solver.

Source: See Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics