EOS algebra is duplicated across call sites; needs a single source of truth

[sbryngelson]

Summary

The stiffened-gas EOS algebra has canonical definitions in src/common/m_variables_conversion.fpp — s_compute_pressure, s_compute_internal_energy, s_compute_speed_of_sound — but several call sites open-code the same expressions instead of calling them. Two have already drifted into real bugs (#1706, #1707). This issue tracks making the algebra DRY, which is prerequisite groundwork for the multi-EOS interface in #1638: a second EOS cannot be introduced safely while the first one has several independent definitions.

Pressure inversion

Canonical, m_variables_conversion.fpp:71:

pres = (energy - dyn_p - pi_inf - qv)/gamma

Open-coded copies:

Site	Expression	Status
simulation/m_pressure_relaxation.fpp:311	(E - dyn_pres - qv_mix - pi_inf)/gamma	correct; operand order differs
simulation/m_data_output.fpp:1558	((E - KE)/(1 - alf) - pi_inf - qv)/gamma	correct; duplicates the bubble branch
simulation/m_data_output.fpp:1621	((E - KE)/(1 - alf) - pi_inf - qv)/gamma	correct; duplicates the bubble branch
simulation/m_bubbles_EL.fpp:374	(E - dynP - pi_inf)/gamma	wrong — omits qv (#1706)

Specific total enthalpy at sound-speed call sites

s_compute_speed_of_sound requires an H that includes qv (it subtracts qv/rho internally). The 27 call sites split into two conventions:

• Solver paths derive it from conserved energy — H_L = (E_L + pres_L)/rho_L — which is correct and cannot drift.
• Four diagnostic sites open-code the closed form instead. post_process/m_data_output.fpp:1289 includes qv and is correct; simulation/m_data_output.fpp:1297,1379,1438 omit it and are wrong (bug: simulation probe sound speed omits qv from enthalpy, disagreeing with post-process #1707).

A fifth site, s_compute_cson_from_pinf in simulation/m_bubbles_EL.fpp:844-846,
reconstructs the energy rather than reading it:

E    = gamma*pinf + pi_inf + 0.5_wp*rhol*dot_product(vel, vel)   ! no qv
H    = (E + pinf)/rhol
cson = sqrt((H - 0.5_wp*dot_product(vel, vel))/gamma)            ! no -qv/rho

This one is correct, and must not be "fixed." It omits qv from the
reconstructed E and from the subtraction, and the two omissions cancel
exactly, leaving c^2 = ((gamma+1)*p + pi_inf)/(rho*gamma) — the same result the
canonical routine produces. Adding qv to only one of the two halves would
introduce a bug where none exists today.

It is listed here because it is a fifth independent copy of the algebra whose
correctness depends on a coincidence between two separate omissions. That is
precisely the fragility this issue is about: the next person to touch it has to
rederive the cancellation to know it is safe.

Why this matters beyond the two bugs

Both defects are the same failure: an expression was copied, the canonical version later gained a term, and the copy did not. The golden suite does not catch either one, because qv = 0 in every case it covers.

Adding Mie-Gruneisen or JWL under #1638 multiplies this. Each open-coded site is a place a second EOS silently does not reach — the selector says mie_gruneisen, and m_bubbles_EL keeps computing a stiffened-gas pressure. Deduplicating first turns "add a backend" into a change with a bounded blast radius.

Scope

Route every open-coded site through the canonical routines. Deliberately not in scope here: moving the routines into a new module, or adding backend dispatch. Those are the next steps in #1638; this issue is only about there being one definition of each operation.

Note that folding the sites in changes floating-point operand order in places (- qv_mix - pi_inf vs - pi_inf - qv), so results move in the last bits. The golden tolerance is 1e-12 for double precision (toolchain/mfc/test/case.py:338), roughly four orders of magnitude above that, so no goldens need regenerating — but the two bug fixes above will legitimately move qv /= 0 results and should be reviewed on their own merits.

[sbryngelson]

Adding a fifth site to the enumeration above, found while checking whether the initial Lagrange bubble pressure in #1706 was really meant to be the ordinary liquid pressure.

s_compute_cson_from_pinf (simulation/m_bubbles_EL.fpp:844-846) reconstructs the energy from the driving pressure instead of reading it:

E    = gamma*pinf + pi_inf + 0.5_wp*rhol*dot_product(vel, vel)
H    = (E + pinf)/rhol
cson = sqrt((H - 0.5_wp*dot_product(vel, vel))/gamma)

My first read was that this is another instance of #1707 — the reconstructed E is missing the qv term that E = gamma*p + pi_inf + qv + rho*KE carries. It is not. The subtraction is missing qv/rho as well, and the two cancel:

canonical:  c^2 = ((gamma*p + pi_inf + qv + rho*K + p)/rho - K - qv/rho)/gamma = ((gamma+1)*p + pi_inf)/(rho*gamma)
here:       c^2 = ((gamma*p + pi_inf +      rho*K + p)/rho - K         )/gamma = ((gamma+1)*p + pi_inf)/(rho*gamma)

Same answer. It is correct and should not be "fixed." Adding qv to the reconstructed E without also restoring the - qv/rho — the obvious patch if someone pattern-matches this against #1707 — would introduce a bug where none exists.

I am flagging it precisely because it is a trap. Its correctness rests on a coincidence between two independent omissions, several lines apart, with nothing saying so. Anyone auditing the sound-speed sites after #1707 lands will find it, see a missing qv, and be one keystroke from breaking it.

That is the argument for this issue in miniature. Counting the open-coded copies across both operations:

Operation	Copies	Correct	Wrong	Correct only by cancellation
Pressure inversion	4	3	1 (#1706)	0
Enthalpy / sound speed	5	1	3 (#1707)	1

Nine hand-written copies of two expressions, four of them wrong and one right for a reason no reader can see locally.

When the deduplication happens, this site should route through the shared routine with the true enthalpy and let the canonical subtraction do its job, replacing the coincidence with something verifiable on sight.