Reacting Roe sound speed: c_sum_Yi_Phi is read uninitialized in HLL and hard-coded to zero in HLLC

[sbryngelson]

c_sum_Yi_Phi reaches s_compute_speed_of_sound_avg as the c_c argument, which uses it as

if (chemistry) then  ! Reacting mixture sound speed
    if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
        c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - H))
    else
        call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)
    end if

It is written in exactly two places: roe_avg (src/simulation/include/inline_riemann.fpp), inside its if (chemistry) block, and src/simulation/m_riemann_solver_hllc.fpp:957, where the hypoelastic variant zeroes it. Two consequences follow.

1. HLL never initializes it

c_sum_Yi_Phi is a private variable of the offload loop in src/simulation/m_riemann_solver_hll.fpp and is written only when avg_state == 1, because roe_avg is the only writer. With chemistry, avg_state = 2 and wave_speeds = 2, line 334 passes an uninitialized value, and abs(c_c) > verysmall reads it.

Fortran does not guarantee short-circuit evaluation of .and., so the read can happen even though avg_state == avg_state_roe is already false. The value is not used — the else branch is taken either way — so this is an uninitialized read rather than a wrong answer. It is the shape that tends to surface as a platform-specific NaN rather than a reproducible failure.

2. The HLLC base path hard-codes zero

src/simulation/m_riemann_solver_hllc.fpp:285-286 passes a literal 0._wp where HLL (line 334) and the HLLC hypoelastic variant (line 1050) pass the real value:

call s_compute_speed_of_sound_avg(pres_R, rho_avg, gamma_avg, pi_inf_R, qv_avg, vel_avg_rms, &
                                  & H_avg, 0._wp, alpha_R, c_avg)

Since abs(c_c) > verysmall gates the branch, HLLC without hypoelasticity always falls back to s_compute_speed_of_sound and never uses the reacting Roe sound speed, while HLL and HLLC-with-hypoelasticity do. So the same case can get a different sound speed depending on which solver and which variant it runs through. If that is deliberate it deserves a comment saying so; if not, it is an inconsistency.

[sbryngelson]

Audit corrections before fixing this.

1. Two HLLC sites hard-code 0._wp, not one: m_riemann_solver_hllc.fpp:288 and :640. The
hypoelastic variant at :1058 passes the real c_sum_Yi_Phi, as does HLL at :337.

2. The only initialization is inside #:if HYPO. c_sum_Yi_Phi = 0.0_wp at
m_riemann_solver_hllc.fpp:961 sits under the #:if HYPO block opened at :887. So the non-
hypoelastic HLLC path never initializes it either; it is only safe because it passes a literal.

3. Lax-Friedrichs carries five dead declarations. m_riemann_solver_lf.fpp:60 declares
Cp_avg, Cv_avg, T_avg, eps, c_sum_Yi_Phi and :108 carries all five in the GPU private() list,
but LF never calls s_compute_average_state, s_compute_chemistry_average_state or
s_compute_speed_of_sound_avg - zero uses outside the declaration and the private clause. They can
be deleted outright.

4. The whole reacting Roe path is untested. Every chemistry test case sets wave_speeds = 1, so
s_compute_average_state is never called in any of them. That means the chemistry branch of the Roe
average, c_sum_Yi_Phi, and the c = sqrt(c_c - (gamma-1)*(vel_sum - H)) branch of
s_compute_speed_of_sound_avg have no coverage at all. That is why HLLC silently disabling it went
unnoticed, and it is the same sibling-cases gap that hid the Roe average bug fixed in #1762.

Adding a chemistry + avg_state = 1 + wave_speeds = 2 case is worth doing before deciding
anything about the 0._wp, since right now nothing would detect a change either way.

[sbryngelson]

Research: the HLL/HLLC asymmetry predates the #1762 refactor on both sides.

• HLLC passed a literal 0._wp in the c_c position before 9974b687 ("Derive the speed of sound
  from the state, not a supplied enthalpy"):
  call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, vel_avg_rms, 0._wp, c_avg, qv_avg)
• HLL passed the real c_sum_Yi_Phi before both 9974b687 and 201a93c7.

So the refactor preserved each call faithfully and did not create the divergence. That removes "did the
EOS work break this" from the question and leaves only whether it was ever intended.

Two things make this hard to answer from the code alone. Nothing documents why HLLC's base path opts
out while HLL and HLLC-with-hypoelasticity opt in. And the branch it gates - c = sqrt(c_c - (gamma - 1)*(vel_sum - H)) - has no test coverage at all: every chemistry case in the suite sets
wave_speeds = 1, so s_compute_average_state is never called in any of them, and the reacting Roe
sound speed is never evaluated on any path.

Recommendation: add the coverage before changing the behaviour. #1762 adds
1D -> Chemistry -> Inert Shocktube -> Reacting Roe Average (avg_state = 1, wave_speeds = 2,
riemann_solver = 1) - HLL rather than HLLC precisely because HLLC's literal zero means it would not
exercise the branch. That establishes a golden for the path that is enabled today. Deciding HLLC's
case afterwards then becomes a change with a baseline to compare against, rather than a change to code
nothing has ever tested.

The uninitialized-read half of this issue is fixed in #1762 (c_sum_Yi_Phi is now zeroed in HLL's
chemistry block); this comment concerns only the 0._wp asymmetry.