bug: simulation probe sound speed omits qv from enthalpy, disagreeing with post-process

[sbryngelson]

Summary

Three s_compute_speed_of_sound call sites in src/simulation/m_data_output.fpp build the specific total enthalpy H without the qv term while still passing the real qv as the final argument. Because the routine subtracts qv/rho internally, the sound speed comes out too low whenever qv /= 0. The equivalent post-process path builds H correctly, so simulation and post-process report different sound speeds for identical state.

The contract

s_compute_speed_of_sound (src/common/m_variables_conversion.fpp:1180) computes, in the default stiffened-gas branch (line 1227):

c = (H - 5.e-1*vel_sum - qv/rho)/gamma

Since E = gamma*p + pi_inf + qv + rho*KE, the specific total enthalpy H = (E + p)/rho includes qv. The internal - qv/rho then cancels it, leaving the correct stiffened-gas result c^2 = ((gamma+1)*p + pi_inf)/(rho*gamma). H must therefore include qv.

Every solver call site satisfies this implicitly by deriving H from the conserved energy, e.g. src/simulation/m_riemann_solver_hll.fpp:255:

H_L = (E_L + pres_L)/rho_L

Those paths are correct and unaffected.

The defect

src/simulation/m_data_output.fpp, lines 1297, 1379, 1438:

call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, &
                              & 0._wp, c, qv)

H is open-coded as ((gamma + 1)*pres + pi_inf)/rho — no qv — but the real qv is still passed. The internal subtraction then has nothing to cancel against:

c^2 = ((gamma+1)*p + pi_inf - qv)/(rho*gamma)      instead of
c^2 = ((gamma+1)*p + pi_inf)/(rho*gamma)

src/post_process/m_data_output.fpp:1289 gets it right:

H = ((gamma + 1._wp)*pres + pi_inf + qv)/rho
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, 0._wp, 0._wp, c, qv)

Impact

Diagnostics only — these sites feed probe output and Ma = maxvel/c. Reported Mach number is too high when qv /= 0, and disagrees with what post-process reports for the same field data. If qv > ((gamma+1)*p + pi_inf) the argument goes negative and the mixture_err floor clamps c to 100*sgm_eps, which would make the reported Mach number meaningless rather than merely wrong.

Dormant when qv = 0 (the default). Triggers for the phase-change and reactive-burn examples that set fluid_pp(i)%qv.

Suggested fix

Either add the + qv term at the three sites, or better, stop open-coding H at call sites entirely — the solver convention H = (E + pres)/rho is both correct and drift-proof. See the tracking issue for deduplication.