Post

Numerical Methods for MHD - Part 3: Riemann Solvers

Posted Updated
By Donghui Son 10 min read

Series Navigation

Table of Contents

  1. Riemann Solvers and Solutions for System Equations

Notation note: Throughout this series, we keep $\mu_0$ explicit. In normalized MHD units, the same formulas reduce to the familiar $\mu_0 = 1$ form after nondimensionalization.

5. Riemann Solvers and Solutions for System Equations

A Riemann problem is a standard building block in finite-volume methods: it models how discontinuous initial states $U_L$ and $U_R$ at a cell interface evolve over time. Because MHD can support up to seven wave families (fast, Alfvén, slow, and contact), solving this problem exactly can be quite involved.

5.1 Exact Riemann Solver

  • Setup

    • We specify initial data with a step discontinuity at $x = 0$:

      \[U(x, 0) = \begin{cases} U_L, & x < 0 \\ U_R, & x > 0 \end{cases}\]

  • Wave Structure in MHD

    • Up to seven waves can emerge:

      • Two fast magnetosonic waves.
      • Two Alfvén (shear) waves.
      • Two slow magnetosonic waves.
      • One contact discontinuity.

  • Procedure

    1. Solve for pressure, velocity, and magnetic field in intermediate states.
    2. Determine wave speeds and states between each wave family.
    3. Apply Rankine–Hugoniot conditions across each wave/discontinuity.

  • Complexity

    • An exact MHD Riemann solver is computationally expensive, requiring iterative solutions for intermediate states.
    • Typically used for benchmark test problems or smaller grids, but rarely for large-scale simulations.

5.2 Approximate Riemann Solvers

Because exact solvers can be impractical in large multidimensional simulations, approximate solvers are more common. They provide a balance between accuracy and efficiency by simplifying the wave structure.

5.2.1 Roe’s Approximate Riemann Solver

  • Idea
    • Linearize the MHD equations around a Roe-averaged state $\tilde{U}$.

    • Compute fluxes by decomposing the discontinuity into characteristic waves with approximate speeds.

      \[F_{i+\frac{1}{2}} = \frac{1}{2} \left[ F(U_L) + F(U_R) \right] - \frac{1}{2} \sum_{k=1}^{m} \tilde{\alpha}_k |\tilde{\lambda}_k| \tilde{r}_k,\]

      where:

      • $\tilde{\lambda}_k$ : eigenvalues at $\tilde{U}$ (Roe-averaged state),
      • $\tilde{r}_k$ : corresponding right eigenvectors,
      • $\tilde{\alpha}_k$ : wave strengths.
  • Pros
    • Can capture contact and shear (Alfvén-type) waves more accurately than simpler solvers (e.g. HLL).
    • Less expensive than the full exact MHD solver.
  • Cons
    • The linearization may yield non-physical states if not carefully handled (possible negative densities/pressures).
    • Implementation is more complex than simpler approximate solvers.

5.2.2 HLL Approximate Riemann Solver

  • Key Idea

    The Harten–Lax–van Leer (HLL) solver replaces the full wave fan with just two effective waves at speeds $S_L$ and $S_R$. It assumes there is a single “star” region between these two bounding waves.

    • Minimal Wave Speed: $S_L$
    • Maximal Wave Speed: $S_R$

    These speeds are usually chosen to enclose all possible characteristic waves emanating from the discontinuity. In MHD, one common approach is:

    \[S_L = \min(u_L - c_{fast,L}, u_R - c_{fast,R}), \quad S_R = \max(u_L + c_{fast,L}, u_R + c_{fast,R}),\]

    where $u_L$ and $u_R$ are the normal flow speeds on the left and right, and $c_{fast,L}$, $c_{fast,R}$ are their respective fast magnetosonic speeds.

  • Flux Formula

    When $S_L \lt 0 \lt S_R$, the solver merges the entire wave fan into a single uniform “star” state between those speeds. The final flux at the interface is:

    \[F_{\text{HLL}} = \frac{S_R F(U_L) - S_L F(U_R) + S_L S_R (U_R - U_L)}{S_R - S_L}.\]

    If $S_L \geq 0$, the entire wave fan is moving to the right, so $F_{\text{HLL}} = F(U_L)$. Conversely, if $S_R \leq 0$, everything moves left, so $F_{\text{HLL}} = F(U_R)$.

  • Advantages:

    • Very simple and robust: one does not need to solve for intermediate states inside the wave fan.
    • Positive-preserving for density and pressure, provided wave speeds are chosen carefully.

  • Limitations:

    • HLL omits all internal waves (e.g., contact, shear, slow waves). This can lead to excessive numerical diffusion in MHD, especially around contact discontinuities or shear flows.

5.2.3 HLLC Approximate Riemann Solver

The HLLC (HLL-Contact) solver extends HLL by recovering the contact discontinuity. While the original HLL solver lumps everything into one star region, HLLC introduces a middle “contact” wave speed $S_*$. This allows the solver to better capture density jumps (contacts) and velocity jumps (shear waves in hydrodynamics).

For classical MHD, however, a direct HLLC extension is more delicate than in pure hydrodynamics. The algebra below should therefore be read mainly as the hydrodynamic template behind contact-restoring solvers; in practical ideal-MHD codes, HLLD is often preferred when the normal magnetic field is nonzero.

  • Wave Speeds
    • Left wave at speed $S_L$,
    • Right wave at speed $S_R$,
    • Contact (middle) wave at speed $S_*$.

    The left and right bounding speeds, $S_L$ and $S_R$, can be estimated similarly to HLL. For hydrodynamics (no magnetic field),

    \[S_L = \min(u_L - c_{s,L}, u_R - c_{s,R}), \quad S_R = \max(u_L + c_{s,L}, u_R + c_{s,R}),\]

    where $c_{s,L}$ and $c_{s,R}$ are sound speeds. For MHD, one can instead use fast magnetosonic speeds. The contact wave speed $S_*$ is derived by imposing pressure continuity across the contact. In simpler hydrodynamics, an algebraic formula yields

    \[S_* = \frac{\rho_R u_R (S_R - u_R) - \rho_L u_L (S_L - u_L) + (p_L - p_R)}{ \rho_R (S_R - u_R) - \rho_L (S_L - u_L) },\]

    though details differ in MHD.

  • Flux Formula

    The interface flux $F_{\text{HLLC}}$ is piecewise-defined, depending on whether the solution is moving left, right, or in the star region:

    \[F_{\text{HLLC}} = \begin{cases} F(U_L), & S_L \geq 0, \\ F^{\ast}_{L}, & S_L \leq 0 \leq S_{\ast}, \\ F^{\ast}_{R}, & S_{\ast} \leq 0 \leq S_R, \\ F(U_R), & S_R \leq 0, \end{cases}\]

    where $F_{L}^{\ast} \text{ and } F_{R}^{\ast}$ are the “star-region” fluxes:

    \[F_{L}^{\ast} = F(U_{L}) + S_{L} (U_{L}^{\ast} - U_{L}), \quad F_{R}^{\ast} = F(U_{R}) + S_{R} (U_{R}^{\ast} - U_{R}).\]

    In those expressions, $U_{L}^{\ast}$ and $U_{R}^{\ast}$ represent the conservative variables in the left and right star states.

    For instance, the contact wave ensures pressure is continuous, $p_{L}^{\ast} = p_{R}^{\ast}$, and the normal velocity is $S_{\ast}$.

  • Intermediate State

    • Density and Velocity:

      \[\rho^* = \rho_L \left( \frac{S_L - u_L}{S_L - S_*} \right), \quad u^* = S_*,\]

      if the star region is coming from the left side. Similar expressions hold on the right side.

    • Pressure:

      The contact discontinuity implies

      \[p^*_L = p^*_R = p^*,\]

      so the star states share a common pressure across the contact wave.

  • Advantages:

    • Captures contact (and shear in pure hydrodynamics) discontinuities better than HLL.
    • Generally less diffusive near contact waves.
  • Limitations:
    • In MHD, HLLC does not fully resolve Alfvén and slow waves.
    • May still introduce moderate diffusion around complex wave interactions in magnetized flows.
    • Certain extreme states can produce unphysical pressure/density if wave speed estimates are inaccurate.

5.2.4 HLLD Approximate Riemann Solver

In ideal MHD, the exact 1D Riemann problem contains seven characteristic waves (two fast, two slow, two Alfvén, and one contact). The HLLD solver does not retain all seven explicitly. Instead, it approximates the fan with a five-wave structure: two outer fast waves, two inner rotational (Alfvén) discontinuities, and one contact discontinuity in the middle. This is the key reason HLLD is both more accurate than HLL/HLLC for MHD and still much cheaper than an exact solver.

Overview
  • Resolved explicitly: left/right fast waves, left/right rotational discontinuities, and the contact wave.
  • Not resolved explicitly: the slow waves are absorbed into the reconstructed intermediate states rather than carried as separate waves.
  • Practical consequence: HLLD sharply improves the treatment of contact and Alfvénic structures compared with HLL, while avoiding a full seven-wave exact solve.
Wave Speeds

For a 1D interface normal to the $x$-direction, define on each side

\[a_i^2 = \frac{\gamma p_i}{\rho_i}, \qquad b_i^2 = \frac{||\mathbf{B}_i||^2}{\mu_0 \rho_i}, \qquad b_{x,i}^2 = \frac{B_{x,i}^2}{\mu_0 \rho_i}.\]

Then a commonly used estimate of the fast magnetosonic speed is

\[c_{f,i}^2 = \frac{1}{2} \left( a_i^2 + b_i^2 + \sqrt{\left(a_i^2 + b_i^2\right)^2 - 4 a_i^2 b_{x,i}^2} \right).\]

The outer bounding speeds are then chosen as

\[S_L = \min(u_{x,L} - c_{f,L},\, u_{x,R} - c_{f,R}), \qquad S_R = \max(u_{x,L} + c_{f,L},\, u_{x,R} + c_{f,R}).\]

HLLD further assumes that the normal velocity and the total pressure

\[p_t = p + \frac{||\mathbf{B}||^2}{2\mu_0}\]

are constant across the central part of the approximate Riemann fan. This gives the contact speed estimate

\[S_M = \frac{ p_{t,R} - p_{t,L} + \rho_L u_{x,L}(S_L - u_{x,L}) - \rho_R u_{x,R}(S_R - u_{x,R}) }{ \rho_L (S_L - u_{x,L}) - \rho_R (S_R - u_{x,R}) }.\]

The first star-state densities are

\[\rho_L^* = \rho_L \frac{S_L - u_{x,L}}{S_L - S_M}, \qquad \rho_R^* = \rho_R \frac{S_R - u_{x,R}}{S_R - S_M}.\]

From these, the rotational-wave speeds are usually written

\[S_L^* = S_M - \frac{|B_x|}{\sqrt{\mu_0 \rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\mu_0 \rho_R^*}},\]

with the 1D constraint that the normal magnetic field $B_x$ is constant across the interface.

Intermediate States

HLLD constructs four intermediate states:

  • $U_L^*$ between the left fast wave and the left rotational discontinuity,
  • $U_L^{**}$ between the left rotational discontinuity and the contact,
  • $U_R^{**}$ between the contact and the right rotational discontinuity,
  • $U_R^*$ between the right rotational discontinuity and the right fast wave.

Conceptually:

  • Across the fast waves, density, normal momentum, tangential velocity, and tangential magnetic field are updated from the HLL-like jump conditions.
  • Across the rotational discontinuities, the density stays fixed while the tangential velocity and tangential magnetic field rotate consistently with Alfvénic jumps.
  • Across the contact, the normal velocity $S_M$ and total pressure $p_t$ remain continuous.

This five-wave construction is the distinctive feature of HLLD: it restores the most important internal MHD structure without attempting a full exact seven-wave decomposition.

Flux Computation

Once the wave ordering

\[S_L \le S_L^* \le S_M \le S_R^* \le S_R\]

and the corresponding states are known, the interface flux is selected piecewise:

\[F_{i+\frac{1}{2}} = \begin{cases} F(U_L), & S_L \ge 0, \\ F_L^*, & S_L \le 0 \le S_L^*, \\ F_L^{**}, & S_L^* \le 0 \le S_M, \\ F_R^{**}, & S_M \le 0 \le S_R^*, \\ F_R^*, & S_R^* \le 0 \le S_R, \\ F(U_R), & S_R \le 0, \end{cases}\]

where

\[F_L^* = F(U_L) + S_L \left(U_L^* - U_L\right), \qquad F_R^* = F(U_R) + S_R \left(U_R^* - U_R\right),\]

and the double-star fluxes are obtained by applying the same Rankine-Hugoniot logic across the rotational waves:

\[F_L^{**} = F_L^* + S_L^* \left(U_L^{**} - U_L^*\right), \qquad F_R^{**} = F_R^* + S_R^* \left(U_R^{**} - U_R^*\right).\]
Summary of Implementation Steps
  1. Estimate the outer wave speeds $S_L$ and $S_R$ from fast magnetosonic bounds.
  2. Compute the contact speed $S_M$ from continuity of normal velocity and total pressure.
  3. Construct the first star states $U_L^\ast$ and $U_R^\ast$ and then the double-star states $U_L^{\ast\ast}$ and $U_R^{\ast\ast}$.
  4. Select the correct flux branch according to the sign of the ordered wave speeds.
Advantages & Limitations

Advantages

  • Better MHD structure: HLLD resolves contact and rotational discontinuities much more sharply than HLL.
  • Good cost/accuracy balance: It is far cheaper than an exact MHD Riemann solver but captures the most important internal MHD waves.
  • Widely used in production codes: It is a standard choice in many modern Godunov MHD solvers.

Limitations

  • More complex than HLL/HLLC: The intermediate-state algebra is substantially longer and easier to implement incorrectly.
  • Slow waves are not explicit branches: Their effect is represented indirectly inside the intermediate states.
  • Can still fail in extreme regimes: Very low densities, low pressures, or highly magnetized states may require positivity fixes or more diffusive fallbacks.

References

Computational MHD, Introduction
mhd riemann-solver roe-solver hll-solver numerical-methods intermediate
This post is licensed under CC BY 4.0 by the author.