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