Numerical Methods for MHD - Part 3: Riemann Solvers

Series Navigation

• Part 1: Introduction and Basic Equations
• Part 2: Numerical Methods and Solution Techniques
• Part 3: Riemann Solvers
• Part 4: Advanced Topics

Table of Contents

1. Riemann Solvers and Solutions for System Equations

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).

• 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 MHD, there can be up to seven characteristic waves (two fast, two Alfvén, two slow, and one contact) separating different states. The HLLD solver extends simpler HLL-type methods to resolve more internal waves, thereby capturing Alfvén and contact discontinuities more accurately than HLL or HLLC.

Overview

• Captures: Fast/slow magnetosonic waves, Alfvén waves, and the contact discontinuity.
• Wave Structure

  • In 1D MHD, the full solution generally contains seven waves:

    • Fast (left-moving)
    • Slow (left-moving)
    • Alfvén (left-moving)
    • Contact
    • Alfvén (right-moving)
    • Slow (right-moving)
    • Fast (right-moving)

  HLLD approximates these seven waves with carefully chosen bounding speeds $\left( S_L, S_R \right)$ and additional intermediate speeds $\left( S_{\ast L}, S_M, S_{\ast R} \right)$ to reconstruct much of the full wave fan.

Wave Speeds

• Fast wave speeds:

  Estimate $S_L$ (left edge) and $S_R$ (right edge) such that they encompass all possible MHD wave speeds emanating from either side of the discontinuity. A common approach:

  \[c_{fast,L} = \sqrt{ \frac{1}{2} \left( c_{s,L}^2 + v_{A,L}^2 + \sqrt{ ( c_{s,L}^2 + v_{A,L}^2 )^2 - 4 c_{s,L}^2 v_{A,L}^2 } \right), }\] \[c_{fast,R} = \sqrt{ \frac{1}{2} \left( c_{s,R}^2 + v_{A,R}^2 + \sqrt{ ( c_{s,R}^2 + v_{A,R}^2 )^2 - 4 c_{s,R}^2 v_{A,R}^2 } \right), }\] \[S_L = \min(u_L - c_{fast,L}, u_R - c_{fast,R}),\] \[S_R = \max(u_L + c_{fast,L}, u_R + c_{fast,R}),\]

  where:

  • $u_L$ and $u_R$: Normal (or x-directed) flow speeds on the left and right.
  • $c_{s,L}$ and $c_{s,R}$: Sound speeds.
  • $v_{A,L}$, $v_{A,R}$: Alfvén speeds.
  • $c_{fast,L}$, $c_{fast,R}$: Fast magnetosonic speeds on each side.

• Intermediate wave speeds

  • Left- and Right-side Alfvén speeds $S_{\ast L} \text{ and } S_{\ast R}$: typically determined from MHD relations for tangential magnetic fields in the star region.
  • Contact wave speed $S_M$: derived by imposing continuity of total pressure and tangential components of the magnetic field across the contact discontinuity.

Detailed Calculation of Wave Speeds

1. Local Alfvén Speeds

   For a 1D configuration with the primary flow in $x$ and tangential fields in $\left( y, z \right)$:

   \[v_{A,L} = \frac{B_{x,L}}{\sqrt{\mu_0 \rho_L}}, \quad v_{A,R} = \frac{B_{x,R}}{\sqrt{\mu_0 \rho_R}}\]

   However, in a more general 1D MHD setting, $B_x$ is not necessarily the entire magnetic magnitude on each side; we also consider tangential fields $B_{y,L}$, $B_{z,L}$, $B_{y,R}$, $B_{z,R}$.

2. Total (Magnetic + Gas) Pressure

   \[p_{t,L} = p_L + \frac{||\mathbf{B}_{L}||^2}{2 \mu_0}, \quad p_{t,R} = p_R + \frac{||\mathbf{B}_{R}||^2}{2 \mu_0}\]

   The difference in total pressures plays a key role in finding the contact wave speed $S_M$.

3. Contact (Middle) Speed

   The so-called “contact” wave is where the normal velocity and total pressure are continuous across that discontinuity. An approximate formula (from Toro’s text or MHD references) gives:

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

   In HLLD, $B_{y}^{\ast} \text{ and } B_{z}^{\ast}$ remain constant across this contact wave in 1D setups.

4. Intermediate densities

   E.g., the left star region:

   \[\rho^{\ast}_{L} = \rho_L \frac{S_L - u_L}{S_L - S_M}\]

   and similarly for the right star region.

Intermediate States

HLLD acknowledges more sub-waves than HLLC, typically constructing two star states on each side ($U_{L}^{\ast}, U_{L}^{\ast \ast}$ on the left; $U_{R}^{\ast}, U_{R}^{\ast \ast}$ on the right), associated with the Alfvén and slow waves. The tangential magnetic field components are often updated by ensuring proper jump conditions across each wave.

1. Tangential Magnetic Field

   In simpler 1D notations:

   \[B_{y,L}^{\ast} = B_{y,L}, \quad B_{z,L}^{\ast} = B_{z,L}\]

   across the contact. More precise formulas might involve partial jumps if the shear/Alfvén waves are resolved separately $\left(S_{\ast L}, S_{\ast R} \right)$.

2. Compute the contact speed $S_M$:

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

3. Iterative or Direct Formulas

   Some HLLD derivations provide direct expressions for these intermediate states; others rely on iterative adjustments to ensure consistent mass/momentum/energy flux across each wave.

Flux Computation

Once the wave speeds $\left( S_L, S_{\ast L}, S_M, S_{\ast R}, S_R \right)$ and the corresponding states $\left( U_{L}^{\ast}, U_{L}^{\ast \ast}, U_{R}^{\ast}, U_{R}^{\ast \ast} \right)$ are found, the interface flux $F_{i+\frac{1}{2}}$ is piecewise-defined:

1. If $S_L \geq 0$

   \[F_{i+\frac{1}{2}} = F(U_L)\]

   (All waves move to the right.)

2. If $S_L \leq 0 \leq S_{\ast L}$ (Left star region)

   \[F_{i+\frac{1}{2}} = F(U_L) + S_L (U^*_L - U_L)\]

3. If $S_{\ast L} \leq 0 \leq S_M$ (Left star $\rightarrow$ contact)

   \[F_{i+\frac{1}{2}} = F^{\ast}_{L} + S_{\ast L} (U^{\ast \ast}_{L} - U^{\ast}_{L})\]

4. If $S_M \leq 0 \leq S_{\ast R}$ (contact $\rightarrow$ right star)

   \[F_{i+\frac{1}{2}} = F^{\ast}_{R} + S_{\ast R} (U^{\ast \ast}_{R} - U^{\ast}_{R})\]

5. If $S_{*R} \leq 0 \leq S_R$ (Right star region)

   \[F_{i+\frac{1}{2}} = F(U_R) + S_R (U^{\ast}_{R} - U_R)\]

6. If $S_R \leq 0$

   \[F_{i+\frac{1}{2}} = F(U_R)\]

   (All waves move to the left.)

Summary of Steps

1. Estimate Wave Speeds:

   • Compute $S_L$, $S_R$, $S_M$, $S_{\ast L}$, $S_{\ast R}$.

2. Compute Intermediate States:

   • Calculate $\rho_{L}^{\ast}$, $\rho_{R}^{\ast}$, $B_{y}^{\ast}$, $B_{z}^{\ast}$.

3. Compute Fluxes:

   • Determine the region where $0$ lies in the wave speeds.
   • Use the appropriate flux formula based on the region.

Summary of Implementation Steps

1. Estimate Wave Speeds

   • Compute $S_L$, $S_R$ from fast magnetosonic speeds.
   • Compute the contact speed $S_M$, plus shear/Alfvén speeds $S_{\ast L}$ and $S_{\ast R}$ if needed.

2. Compute Intermediate States

   • $\rho_{L}^{\ast}$, $\rho_{R}^{\ast}$, partial velocities, tangential fields, total pressures, etc.
   • Possibly sub-star states $\left( U_{L}^{\ast}, U_{L}^{\ast \ast}, U_{R}^{\ast}, U_{R}^{\ast \ast} \right)$ to account for Alfvén wave splitting.

3. Piecewise Flux

   • Evaluate which wave region the solution belongs to and calculate $F_{i+\frac{1}{2}}$ accordingly.

Advantages & Limitations

Advantages

• Wave Resolution: HLLD can resolve fast, slow, Alfvén, and contact discontinuities more sharply than HLL or HLLC.
• Efficiency: Much faster than a full exact MHD solver, while retaining key wave physics.
• Wide applicability: Suitable for problems with strong shocks, discontinuities, or shear flows in magnetized plasmas.

Limitations

• Complex Implementation: More elaborate than HLL or HLLC, requiring correct wave speed estimates and sub-star state calculations.
• Computational Cost: Higher than simpler solvers (e.g., HLL), though typically still feasible in large-scale MHD codes (e.g., Athena++, PLUTO).
• Approximation: Like all approximate solvers, it can introduce some numerical diffusion or mild instabilities in extreme states (e.g., highly relativistic or extremely low plasma-beta regimes).