Riemann Solver

What Is a Riemann Solver? The Mathematical Structure of Conservation Laws and Wave Propagation

In computational fluid dynamics (CFD) and computational magnetohydrodynamics (CMHD), one of the central tools for understanding the dynamical behavior of physical systems is the Riemann solver. It is essential for numerically solving conservation laws written as systems of hyperbolic partial differential equations. In differential form, a conservation law is typically written as

\[\frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{U}) = 0\]

where $\mathbf{U}$ is the vector of conserved variables, such as mass density $\rho$, momentum density $\rho \mathbf{v}$, total energy density $E$, and magnetic field $\mathbf{B}$, and $\mathbf{F}(\mathbf{U})$ is the corresponding flux.

For clarity, this post focuses primarily on the one-dimensional face-normal Riemann problem, using the ideal-gas Euler equations and ideal MHD as representative examples. Unless stated otherwise, the MHD wave-speed formulas are written in normalized units with $\mu_0 = 1$ or with the factor $4\pi$ absorbed into the nondimensionalization (LeVeque, 2002; Toro, 2009).

Weak Solutions and the Mathematical Structure of Hyperbolic Systems

1. Weak Solutions

An important feature of nonlinear hyperbolic systems is that even smooth initial data can develop infinite gradients in finite time, producing discontinuities, i.e. shock waves. To treat such discontinuous solutions rigorously, one introduces the notion of a weak solution, which satisfies the conservation law in integral form rather than pointwise differential form. For an arbitrary control volume $\Omega$ and time interval $[t_1, t_2]$, this takes the form

\[\int_{\Omega} \mathbf{U}(x, t_2)\, d\Omega - \int_{\Omega} \mathbf{U}(x, t_1)\, d\Omega + \int_{t_1}^{t_2} \oint_{\partial \Omega} \mathbf{F}(\mathbf{U}) \cdot \mathbf{n}\, dS\, dt = 0\]

Weak solutions make it possible to define the solution even when shock surfaces and other discontinuities are present.

2. Hyperbolicity and Wave Propagation

Saying that a system is hyperbolic means that information propagates as waves with finite speed. In one spatial dimension, if the Jacobian matrix

\[A(\mathbf{U}) = \frac{\partial \mathbf{F}}{\partial \mathbf{U}}\]

has real eigenvalues $\lambda_k$ and a complete set of eigenvectors $r_k$ together with left eigenvectors $l_k$, then the system is hyperbolic.

• Eigenvalues $\lambda_k$ represent the characteristic speeds of the wave families.
• Right eigenvectors $r_k$ describe the structure of the corresponding waves, i.e. how the conserved variables vary along each wave family.
• Left eigenvectors $l_k$ are used to measure wave strengths, with $l_i r_j = \delta_{ij}$.

3. Characteristic Decomposition

Using the factorization

\[A = R \Lambda L,\]

where $\Lambda = \operatorname{diag}(\lambda_k)$, the columns of $R$ are the right eigenvectors, and the rows of $L$ are the left eigenvectors, one can introduce characteristic variables

\[d\mathbf{W} = L\, d\mathbf{U}.\]

Locally, the system can then be viewed as $m$ scalar advection equations,

\[\frac{\partial w_k}{\partial t} + \lambda_k \frac{\partial w_k}{\partial x} = \text{(nonlinear coupling terms)},\]

which makes the wave-propagation nature of the dynamics explicit.

4. Riemann Invariants

For some systems, there exist combinations of variables that remain constant along particular characteristic curves $dx/dt = \lambda_k$. These are the Riemann invariants.

• Example: the one-dimensional isentropic Euler equations with $p = K \rho^\gamma$

  The eigenvalues are

  \[\lambda_{1,2} = u \mp c_s,\]

  where $c_s = \sqrt{dp/d\rho}$ is the sound speed. The corresponding Riemann invariants are

  \[R_\mp = u \mp \int \frac{c_s(\rho)}{\rho}\, d\rho = u \mp \frac{2c_s}{\gamma - 1}, \qquad \gamma \ne 1.\]

  These remain constant along the corresponding characteristic curves $\lambda_{1,2}$ and greatly simplify the construction of rarefaction-wave solutions (LeVeque, 2002; Toro, 2009).

The Riemann Problem: Modeling an Initial Discontinuity

The Riemann problem is one of the most fundamental initial-value problems for conservation laws. In one dimension, it consists of solving

\[\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}(\mathbf{U})}{\partial x} = 0\]

subject to discontinuous step data

\[\mathbf{U}(x,0) = \begin{cases} \mathbf{U}_L, & x < 0, \\ \mathbf{U}_R, & x > 0. \end{cases}\]

Here $\mathbf{U}_L$ and $\mathbf{U}_R$ are constant states. The solution often has a self-similar form, meaning that

\[\mathbf{U}(x,t) = \widetilde{\mathbf{U}}(\xi), \qquad \xi = \frac{x}{t}.\]

The Structure of the Exact Riemann Solution

The solution $\mathbf{U}(x,t) = \widetilde{\mathbf{U}}(x/t)$ forms a wave fan in the $x$-$t$ plane, consisting of elementary waves emanating from the origin. An $x$-$t$ diagram is often the clearest way to visualize this structure.

1. Elementary Wave Types

1. Rarefaction wave: A continuous region in which the state varies smoothly with $\xi = x/t$, satisfying $\lambda_k(\mathbf{U}) = \xi$ and $d\mathbf{U}/d\xi \propto r_k$.

2. Shock wave: A discontinuous jump. If the shock moves with speed $s$, then the left and right states $\mathbf{U}_L$ and $\mathbf{U}_R$ must satisfy the Rankine-Hugoniot condition

   \[\mathbf{F}(\mathbf{U}_R) - \mathbf{F}(\mathbf{U}_L) = s(\mathbf{U}_R - \mathbf{U}_L).\]

   This condition alone does not guarantee uniqueness, because nonphysical expansion shocks may also satisfy it. One therefore imposes an entropy condition. For a $k$-shock in a strictly hyperbolic system, the Lax entropy condition requires

   \[\lambda_k(\mathbf{U}_L) > s > \lambda_k(\mathbf{U}_R).\]

   When neighboring wave families exist, one also requires

   \[\lambda_{k-1}(\mathbf{U}_L) < s < \lambda_{k+1}(\mathbf{U}_R),\]

   so that no additional wave family is artificially compressed across the shock (LeVeque, 2002; Toro, 2009). This selection is closely related to the vanishing-viscosity limit and, physically, to the second law of thermodynamics.

3. Contact discontinuity: For the one-dimensional Euler equations, the pressure $p$ and normal velocity $u$ remain continuous across the contact, while the density $\rho$ and entropy may jump. The discontinuity therefore moves with speed $\xi = u$. In ideal MHD, one must distinguish this from a tangential discontinuity, which can arise when $B_n = 0$; the two should not be conflated (LeVeque, 2002; Toro, 2009).

2. Phase-Plane Analysis for the Euler Equations ($p$-$u$ Plane)

For the ideal-gas Euler equations, provided no vacuum forms, the exact Riemann solution can be constructed in the $p$-$u$ phase plane by finding the intersection $(p^\ast, u^\ast)$ of the left 1-wave curve $u = u_L^\ast(p)$ and the right 3-wave curve $u = u_R^\ast(p)$. The two star states on either side of the contact discontinuity share the same pressure and normal velocity, namely $p^\ast$ and $u^\ast$ (Toro, 2009).

This construction becomes especially clear when shown in a $p$-$u$ phase-plane diagram.

In practice, the intersection is obtained by solving a nonlinear scalar equation such as

\[f(p) = u_L^\ast(p) - u_R^\ast(p) = 0\]

with Newton-Raphson or a similar iteration. If the initial states are separated strongly enough, a vacuum state can form, in which case a separate vacuum branch must be considered (Toro, 2009).

3. Ideal-MHD Wave Structure and Eigenvalues

The face-normal one-dimensional Riemann problem in ideal MHD typically has seven wave families: two fast magnetosonic waves, two Alfvén waves, two slow magnetosonic waves, and one entropy/contact wave. Even in a three-dimensional simulation, the local Riemann problem at a cell face is still usually interpreted through this seven-wave fan.

By contrast, the so-called eight-wave picture does not arise simply because the simulation is three-dimensional. Rather, it appears when the governing equations are augmented so that $\nabla \cdot \mathbf{B}$ errors are carried as an additional divergence mode, as in Powell-type formulations or GLM/divergence-cleaning approaches (Powell et al., 1999; Dedner et al., 2002).

The classical eight-wave picture is associated with Powell-type formulations. By contrast, Dedner-style GLM-MHD introduces an auxiliary scalar cleaning field and typically leads to a nine-wave hyperbolic structure.

In normalized units, the fast and slow magnetosonic speeds are

\[c_{f,s}^2 = \frac{1}{2} \left[ (c_s^2 + v_A^2) \pm \sqrt{(c_s^2 + v_A^2)^2 - 4 c_s^2 v_{Ax}^2} \right],\]

where

\[v_A^2 = \frac{|\mathbf{B}|^2}{\rho}, \qquad v_{Ax}^2 = \frac{B_x^2}{\rho}.\]

Because the eigenstructure is intricate and can become degenerate, implementing an exact MHD Riemann solver is both difficult and computationally expensive (Miyoshi & Kusano, 2005).

The Role of the Riemann Solver in the Finite-Volume Method

In the finite-volume method, the cell average $\mathbf{U}_i^n$ is updated according to

\[\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x} \left( \mathbf{F}_{i+1/2} - \mathbf{F}_{i-1/2} \right).\]

The central role of the Riemann solver is to compute the numerical flux $\mathbf{F}_{i+1/2}$ at each cell interface $x_{i+1/2}$. This is usually done by solving, exactly or approximately, the local Riemann problem defined by the left and right interface states $\mathbf{U}_{i+1/2,L}$ and $\mathbf{U}_{i+1/2,R}$. These states may be taken directly from the neighboring cell averages in a first-order scheme, or reconstructed more accurately from within each cell using methods such as MUSCL or WENO.

For an explicit Godunov-type finite-volume scheme, one typically requires a CFL condition

\[\sigma = \max(|\lambda_k|)\frac{\Delta t}{\Delta x} \le \sigma_{\max}.\]

For one-dimensional first-order monotone schemes, one often takes $\sigma_{\max} \le 1$, but the actual admissible value may be smaller depending on the reconstruction, time integrator, spatial dimension, and source-term treatment (LeVeque, 2002; Toro, 2009).

In multiple dimensions, the standard approach is to solve a local one-dimensional Riemann problem in the face-normal direction at each interface. The resulting normal fluxes can then be used either in an unsplit update or within an operator-splitting framework.

Approximate Riemann Solvers: Ideas and Main Types

Because exact solvers are often too expensive, a wide variety of approximate Riemann solvers have been developed.

1. Godunov’s Method

Godunov’s method defines the interface flux by solving the exact or approximate Riemann problem at the interface and setting

\[\mathbf{F}_{i+1/2} = \mathbf{F}(\widetilde{\mathbf{U}}(0)).\]

It is first-order accurate and serves as the conceptual foundation for many other solvers.

2. Roe Solver

• Principle: Use a locally linearized Roe matrix $\widetilde{A}$ satisfying, among other properties,

  \[\mathbf{F}_R - \mathbf{F}_L = \widetilde{A}(\mathbf{U}_R - \mathbf{U}_L).\]

• Flux formula:

  \[\mathbf{F}_{\text{Roe}} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{1}{2}\sum_k |\widetilde{\lambda}_k|\, \widetilde{\alpha}_k \widetilde{r}_k.\]

  This can be viewed as a central flux plus a characteristic-based numerical viscosity term.

• Characteristics: Roe’s Property U concerns a linearization that preserves the essential upwind structure of the exact Riemann solution; it is not a condition that guarantees positivity of density or pressure (Roe, 1981). The standard Roe solver resolves isolated contact discontinuities very well, but it requires an entropy fix in transonic rarefaction problems. In sufficiently strong rarefactions, density or pressure may fail to remain positive, so positivity-preserving modifications or HLL-type fallbacks are often used in practice. Roe averages also become more delicate for complicated equations of state and for MHD.

3. HLL-Family Solvers (HLL, HLLC, HLLD)

• Principle: Approximate the entire Riemann fan by using only the fastest left- and right-going wave speeds $S_L$ and $S_R$.

• HLL flux:

  \[\mathbf{F}_{\text{HLL}} = \begin{cases} \mathbf{F}_L, & \text{if } S_L \ge 0, \\ \dfrac{S_R \mathbf{F}_L - S_L \mathbf{F}_R + S_L S_R(\mathbf{U}_R - \mathbf{U}_L)}{S_R - S_L}, & \text{if } S_L < 0 < S_R, \\ \mathbf{F}_R, & \text{if } S_R \le 0. \end{cases}\]

• HLLC: Improves HLL by restoring the contact discontinuity through the additional middle speed $S_\ast$, which requires star-state definitions $\mathbf{U}_L^\ast$ and $\mathbf{U}_R^\ast$.

• HLLD: For ideal MHD, HLLD restores additional internal structure, including Alfvénic discontinuities, through a multi-state construction (Miyoshi & Kusano, 2005).

• Characteristics: Very robust overall and comparatively straightforward to implement. Plain HLL is more diffusive, while HLLC and HLLD reduce that diffusion significantly.

4. Flux Vector Splitting (FVS: Steger-Warming, van Leer)

• Principle: Split the Jacobian $A$ or flux $\mathbf{F}$ into positive and negative parts according to the signs of eigenvalues or Mach-number-based factors, then upwind the two parts separately:

  \[\mathbf{F}_{\text{FVS}} = \mathbf{F}^+(\mathbf{U}_L) + \mathbf{F}^-(\mathbf{U}_R).\]

• Characteristics: Very robust and easy to implement, but often too diffusive, especially for contact discontinuities and shear layers.

5. AUSM-Family Solvers (AUSM, AUSM+, AUSMDV)

• Principle: Split the flux into convective and pressure parts and upwind them differently.

• Characteristics: A compromise between the robustness of FVS and the sharpness of flux-difference splitting methods such as Roe. Improved variants are especially good at handling contact discontinuities.

Additional Considerations in MHD

In MHD simulations, several extra issues become especially important:

• One must clearly distinguish between the seven-wave structure of the face-normal ideal-MHD Riemann problem and the additional divergence mode introduced in Powell- or GLM-type formulations.

• One must also maintain the magnetic divergence constraint

  \[\nabla \cdot \mathbf{B} = 0.\]

  Common approaches include:

  • Powell’s 8-wave formulation: adds source terms so that divergence errors are transported as an additional wave mode (Powell et al., 1999).
  • Hyperbolic divergence cleaning: introduces an auxiliary scalar field that propagates and damps divergence errors (Dedner et al., 2002).
  • Constrained transport (CT): stores magnetic fields on faces and electric fields on edges so that Faraday’s law preserves the discrete divergence to roundoff or machine precision on a staggered grid.

Numerical-Analysis Viewpoint

• Order of accuracy: Local truncation error and global error are controlled by the spatial reconstruction and time integration. High-order schemes typically combine reconstruction methods such as MUSCL or WENO with time integrators such as SSP Runge-Kutta methods.
• Stability analysis: One studies both linear stability, for example through von Neumann analysis and CFL restrictions, and nonlinear stability, through ideas such as TVD/TVB behavior and slope limiters.

Summary Comparison of Solver Families

Property	Exact Solver	Roe Solver	HLL Family (HLL/HLLC/HLLD)	FVS Family	AUSM Family
Internal wave structure represented	Resolves the full wave fan	Decomposes all locally linearized propagation modes	HLL keeps only outer waves; HLLC/HLLD recover part of the internal structure	Reflected indirectly through upwind flux splitting	Captures major wave behavior through convective/pressure splitting
Resolution of contact/shear structures	Exact	Excellent	Diffusive in HLL, much improved in HLLC/HLLD	Low	Strong
Positivity / robustness	Exact solution but costly and intricate to implement	Standard form does not guarantee positivity and may need fixes	Generally robust; HLL-type schemes are often favorable for positivity	Very robust but highly diffusive	Generally robust
Applicability to MHD	Possible but very difficult	Possible, but Roe averaging and entropy treatment are delicate	HLLD is especially popular for ideal MHD	Usable but rather diffusive	Baseline forms target Euler flow; MHD requires specialized extensions
Computational cost	Very high	Moderate to high	Low to moderate	Low	Low to moderate

This table compares the intrinsic tendencies of the Riemann solvers themselves. The global order of accuracy for smooth solutions is determined mainly by the reconstruction and time integration rather than by the flux function alone.

Conclusion and Further Discussion

The Riemann solver is the mathematical and numerical bridge between the analytical structure of hyperbolic conservation laws and practical computational simulation. Its theory draws together characteristic waves, discontinuous solutions, exact self-similar wave fans, approximate upwind constructions, and the stability requirements of finite-volume methods.

Each solver family embodies a different trade-off among accuracy, robustness, and computational cost. The best choice depends on the problem at hand: whether contact discontinuities must be sharply resolved, whether strong MHD wave interactions are important, whether positivity must be preserved aggressively, and how much computational expense is acceptable.

In real applications, source terms are often present as well. This leads naturally to more advanced numerical ideas such as well-balanced schemes, which preserve important steady states, and entropy-stable schemes, which aim to control nonlinear stability over long integrations. These topics point beyond classical Riemann solvers toward the broader modern theory of robust shock-capturing methods.

References

• R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
• E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, 2009.
• E. Godlewski and P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer, 1996.
• P. L. Roe, “Approximate Riemann solvers, parameter vectors, and difference schemes”, Journal of Computational Physics 43(2), 357-372, 1981.
• K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw, “A Solution-Adaptive Upwind Scheme for Ideal Magnetohydrodynamics”, Journal of Computational Physics 154(2), 284-309, 1999.
• A. Dedner, F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer, and M. Wesenberg, “Hyperbolic Divergence Cleaning for the MHD Equations”, Journal of Computational Physics 175(2), 645-673, 2002.
• T. Miyoshi and K. Kusano, “A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics”, Journal of Computational Physics 208(1), 315-344, 2005.