# Error in HLL and HLLC Riemann Solvers

**URL:** <https://pyfr.discourse.group/t/error-in-hll-and-hllc-riemann-solvers/1211>\
**Category:** General\
**Created:** [5 February 2025 22:20 UTC](https://pyfr.discourse.group/t/error-in-hll-and-hllc-riemann-solvers/1211 "2025-02-05T22:20:07Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Hasti](https://avatars.discourse-cdn.com/v4/letter/h/41988e/32.png) [@Hasti](https://pyfr.discourse.group/u/Hasti)\
**Post date:** [5 February 2025 22:20 UTC](https://pyfr.discourse.group/t/error-in-hll-and-hllc-riemann-solvers/1211/1 "2025-02-05T22:20:07Z")

</div>

Hello,

Recently I was looking at the Riemann solvers in PyFR and I believe I’ve spotted an error in the HLL and HLLC solvers. The errors are regarding the roe averaged speed of sound calculations.

In the HLL solver:  
pyfr/solvers/euler/kernels/rsolvers/hll.mako:

```
// Get the normal left and right velocities
fpdtype_t nvl = ${pyfr.dot('n[{i}]', 'vl[{i}]', i=ndims)};
fpdtype_t nvr = ${pyfr.dot('n[{i}]', 'vr[{i}]', i=ndims)};

// Compute the Roe-averaged velocity
fpdtype_t nv = (sqrt(ul[0])*nvl + sqrt(ur[0])*nvr)
             / (sqrt(ul[0]) + sqrt(ur[0]));
...
// Roe average sound speed
fpdtype_t a = sqrt(${c['gamma'] - 1}*(H - 0.5*nv*nv));

```

‘nv’ seems to be the normal component of velocity.

And in the HLLC solver:  
pyfr/solvers/euler/kernels/rsolvers/hllc.mako:

```
${pyfr.expand('inviscid_flux_1d', 'ul', 'fl', 'pl', 'vl')};
${pyfr.expand('inviscid_flux_1d', 'ur', 'fr', 'pr', 'vr')};
...
// Roe average sound speed
fpdtype_t u = (sqrt(ul[0])*vl[0] + sqrt(ur[0])*vr[0]) /
              (sqrt(ul[0]) + sqrt(ur[0]));
fpdtype_t a = sqrt(${c['gamma'] - 1}*(H - 0.5*u*u));

```

‘u’ seems to be the x component of velocity.

As far as I understand, the total magnitude of the velocity vector should be used when calculating the Roe average speed of sound, not just one component of velocity. Looking at the Roe solver:

pyfr/solvers/euler/kernels/rsolvers/roe.mako:

```
${pyfr.expand('inviscid_flux_1d', 'ul', 'fl', 'pl', 'vl')};
${pyfr.expand('inviscid_flux_1d', 'ur', 'fr', 'pr', 'vr')};
...
% for i in range(ndims):
    va[${i}] = (vl[${i}]*sqrt(ul[0]) + vr[${i}]*sqrt(ur[0])) * inv_rar;
% endfor

fpdtype_t qq = ${pyfr.dot('va[{i}]', i=ndims)};
fpdtype_t a = sqrt(${c['gamma'] - 1}*(ha - 0.5*qq));

```

The magnitude of the velocity vector is used.

In Toro\* pg 328, when discussing HLL and HLLC solvers, it’s mentioned that Davis\*\*\* and Einfeldt\*\*\* proposed using the Roe average eigenvalues for the left and right waves and the following equation is given:

a = sqrt( (gamma-1) \* (H - 0.5\*u\*u) )

However looking at pg 356 where the Roe solver is covered, the following equation is used:

a = sqrt( (gamma-1) \* (H - 0.5\*V\*V) )  
where V\*V = u\*u + v\*v + w\*w.

Looking further, Einfeldt’s\*\* paper appears to only consider a 1D Euler system, and Davis\*\*\* never clearly states an equation of the speed of sound.

Looking at other open source CFD solvers, namely SU2, the magnitude of the velocity vector is used, not just a velocity component. (line 164 [SU2/SU2\_CFD/src/numerics/flow/convection/hllc.cpp at 0e1495c7ff40afefc7d8613d110115fca40d6782 · su2code/SU2 · GitHub](https://github.com/su2code/SU2/blob/0e1495c7ff40afefc7d8613d110115fca40d6782/SU2_CFD/src/numerics/flow/convection/hllc.cpp#L134)).

From a physical perspective, using the magnitude of the velocity seems to be the most appropriate method in my eyes and I struggle to find a physical justification for using only one component of the velocity.

Many thanks,  
Hasti.

\* Riemann Solvers and Numerical Methods for Fluid Dynamics, Toro  
\*\* On Godunov-Type Methods for Gas Dynamics, Einfeldt  
\*\*\* Simplified Second-Order Godunov-Type Methods, Davis

---

<div class="post-metadata">

**Author:** ![fdw](https://avatars.discourse-cdn.com/v4/letter/f/a587f6/32.png) [@fdw](https://pyfr.discourse.group/u/fdw)\
**Post date:** [5 February 2025 23:35 UTC](https://pyfr.discourse.group/t/error-in-hll-and-hllc-riemann-solvers/1211/2 "2025-02-05T23:35:24Z")

</div>

Both implementations use the magnitude of the velocity vector. Note that the HLLC routine first performs a coordinate transform so \mathbf{n} = (1, 0, 0). This is why it only needs to consider the x component.

Regards, Freddie.