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

Activity

1. sbryngelson commented on Aug 8, 2026

   

   sbryngelson

   on Aug 8, 2026

   Last edited by sbryngelson

   MemberAuthor

   More actions

   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.
2. sbryngelson commented on Aug 26, 2026

   

   sbryngelson

   on Aug 26, 2026

   MemberAuthor

   More actions

   Audit of where the EOS actually lives, surveyed by parameter array (gammas, pi_infs, qvs, gs_min, ps_inf, cvs, qvps). Searching for a guessed expression shape missed whole subsystems twice, so this starts from the parameters and assigns every file that consumes one.

   It is ten expression families, not two: sound speed, mixture coefficients, energy, pressure inversion, enthalpy, bulk modulus, temperature inversion, entropy, caloric relations, and the Tait-named forms (which turn out to be stiffened gas under different names, not a separate EOS).

   #1762 closes three of them — mixture coefficients, energy, sound speed — across the five Riemann solvers, CBC, the viscous fluxes and the acoustic source. Zero sound-speed computations remain outside the operator.

   Still open: 13 files, 96 parameter references.

   file	refs	family
   m_phase_change	39	entropy, temperature, caloric; sole consumer of cvs/qvps
   m_igr	20	mixture coefficients (10 copies), energy
   m_pressure_relaxation	14	pressure inversion, mixture coefficients, isentrope
   m_rhs, m_hypoelastic	6	bulk modulus with a shear term
   m_bubbles_EE, m_qbmm	6	Tait-form isentrope
   pre_process (m_assign_variables, m_data_output)	5	initial conditions
   m_ibm	2	energy from an image point
   post_process (m_derived_variables)	2	derived fields
   m_sim_helpers	1	energy
   m_reactive_burn, m_start_up	2	temperature inversion

   m_phase_change alone is a comparable body of EOS algebra to everything #1762 changed.

   One constraint found along the way that is easy to miss: adv does not mean the same thing in every model. Under bubbles_euler with num_fluids == 1, eqn_idx%alf == eqn_idx%adv%end, so the sole advection slot is the void fraction and the mixture rule does not apply — the coefficients are the liquid's. Any operator taking adv has to go through something that knows this.
3. sbryngelson commented on Aug 28, 2026

   

   sbryngelson

   on Aug 28, 2026

   MemberAuthor

   More actions

   Audit against the current fix/speed-of-sound-from-state branch: every site listed here is now routed
   through the canonical operators.

   site in this issue	current state
   m_pressure_relaxation open-coded pressure	f_pressure(...) at :220
   m_data_output bubble-branch pressure (×2)	s_compute_pressure(...), six call sites
   m_bubbles_EL:374 — omitted qv (#1706)	f_pressure(..., gamma, pi_inf, qv) — the qv term is present
   m_data_output enthalpy sites (#1707)	no open-coded H = ... remains
   s_compute_cson_from_pinf — must not change	intact; still omits qv from both the reconstructed E and the subtraction, so the cancellation stands

   That last one is worth calling out: m_bubbles_EL.fpp was modified on this branch, so the routine was
   edited around deliberately rather than left alone by luck.

   Since this issue was filed, the same treatment has been extended past its original list:

   • m_igr's CFL routine was the last place writing the stiffened-gas pressure and sound speed inline;
     it now uses f_pressure and f_bulk_modulus.
   • The QBMM Keller-Miksis liquid sound speed was n_tait*(pres + B_tait), which is f_bulk_modulus
     written out.
   • The HLLC six-equation star pressure applied the Tait isentrope inline; it is now
     f_pressure_on_isentrope.
   • The CBC mixture-coefficient time derivative had no counterpart at all; it is now
     s_compute_mixture_coefficients_dt, mirroring s_compute_mixture_coefficients branch for branch.
   • The stiffened-gas thermal law p + B = (n-1) cv rho T was written out separately in
     m_phase_change and m_reactive_burn; both now call f_sg_thermal.
   • The post-process gas internal energy dropped the pi_inf and qv terms; it now uses
     f_phase_internal_energy (post_process gas internal energy diagnostic drops pi_inf and qv, and hard-codes fluid 2 #1781).

   Your tolerance analysis was right and mine was wrong. This issue predicted that folding the sites in
   moves results only in the last bits, and that the 1e-12 golden tolerance is four orders of magnitude
   above that, so nothing needs regenerating. I predicted golden movement twice while doing this work - once
   for the IGR substitution and once for QBMM - and both times every golden held. Reading this issue's
   analysis first would have saved both wrong calls.

   Still deliberately outside the scope: m_phase_change's caloric relations and Newton solve,
   m_pressure_relaxation's isentrope inversion, and a few single-fluid gammas(1) diagnostic
   shortcuts. Those need a thermodynamically complete EOS interface rather than deduplication.
4. sbryngelson commented on Sep 2, 2026

   

   sbryngelson

   on Sep 2, 2026

   MemberAuthor

   More actions

   Fixed by #1762, merged as c9c8b0f. This issue was listed in that PR's Closes line but GitHub only auto-closes the first reference in a comma-separated list, so it stayed open; closing manually.
5. sbryngelson

   closed this as completedon Sep 2, 2026