Lagrange bubble initialization silently overwrites the user-supplied global pref

[sbryngelson]

Summary

s_add_bubbles assigns the module-level global pref — a documented, user-settable case parameter — from the initial gas pressure of whichever Lagrangian bubble happens to be added last. The user's configured value is silently discarded.

Location

src/simulation/m_bubbles_EL.fpp:381-384

gas_p(bub_id, 1) = pliq + 2._wp*(1._wp/Web)/bub_R0(bub_id)
if (.not. f_approx_equal((1._wp/Web), 0._wp)) then
    pref = gas_p(bub_id, 1)
end if

There is no local pref in m_bubbles_EL.fpp; this writes the global declared alongside rhoref and defaulted in m_global_parameters_common.fpp:378.

Why this looks wrong

• pref is a user input. It is registered at toolchain/mfc/params/definitions.py:739, described as "Reference pressure", documented in case.md (line 137, and line 876 as "Reference pressure for bubble models"), and required by model_eqns = 4 (definitions.py:405). Overwriting it discards a value the user set.

• It is per-bubble, in a loop. s_add_bubbles is called once per bubble, so the surviving value is whatever the last-added bubble produced. It has no defined meaning — with one bubble it is that bubble's pressure, with a cloud it is an arbitrary member's.

• It is read by the EOS. s_compute_pressure uses it in the four-equation Tait branch (m_variables_conversion.fpp:105):

  pres = (pref + pi_inf)*(energy/(rhoref*(1 - alf)))**(1/gamma + 1) - pi_inf

  It is also pushed to the device in the GPU_UPDATE list at m_global_parameters.fpp:916, so the overwritten value reaches GPU kernels.

• It is conditional on surface tension. The assignment only happens when 1/Web is not approximately zero, so whether a user's pref survives depends on whether surface tension is active — a coupling with no evident justification.

Impact

Unclear, and that is the problem. Whether it currently causes wrong results depends on whether any supported configuration combines Lagrangian bubbles, nonzero surface tension, and a later consumer of pref. The four-equation model is the obvious candidate. Even where nothing reads it, a user-supplied parameter being silently replaced is a trap.

How this surfaced

While reviewing #1709. That PR corrects pliq by including the qv term, which changes gas_p, which changes whatever pref gets overwritten with. The fix is right on its own terms, but it made this assignment visible.

Suggested direction

If the intent is a bubble-model reference pressure derived from initial conditions, it should be its own module variable with a name that says so, not the user-facing pref. If the intent is to honour a user-supplied pref, this assignment should be deleted. Either way it should not depend on whether surface tension is enabled.

Flagging rather than fixing, since the intent is not recoverable from the code and the right resolution depends on what the EL model is supposed to do here.

[sbryngelson]

Traced every read and write of pref. The picture is sharper than "intent unrecoverable", and it points at a name mix-up rather than a design decision.

Every use of the global

Writes

Site	Condition
m_global_parameters_common.fpp:378	default, dflt_real
pre_process/m_global_parameters.fpp:520	bubbles_euler .and. polytropic -> 1._wp
post_process/m_global_parameters.fpp:419	same
simulation/m_bubbles_EL.fpp:381	bubbles_lagrange + surface tension -> last bubble's gas_p

Reads

Site	Guard
m_variables_conversion.fpp:105	model_eqns == 4 (Tait inversion)
pre_process/m_assign_variables.fpp:409	else of if (model_eqns /= model_eqns_4eq) — so model_eqns == 4
pre_process/m_icpp_patches.fpp:646	explicitly model_eqns == model_eqns_4eq
2dHardcodedIC.fpp:115,121, 3dHardcodedIC.fpp:158,164	hardcoded-IC patches, hydrostatic ambient

Consequence

Every EOS-side consumer of pref is gated on model_eqns = 4, and model_eqns = 4 is never exercised: the test suite only ever sets model_eqns to 2 or 3, and no example or benchmark sets 4 either. No case anywhere sets pref.

So the only live read is the hardcoded-IC hydrostatic use in pre-process, which is case-specific by design and happens before the simulation-side write ever occurs.

That makes m_bubbles_EL.fpp:381 a write with no reachable reader in any exercised configuration. It is currently harmless — but it is harmless by accident, not by construction, and it is still GPU_UPDATEd to device.

The likely explanation

The Euler-Lagrange model has its own reference set, declared separately from pref/rhoref:

real(wp) :: p0ref, rho0ref, T0ref, ss, pv, vd, ...     ! m_global_parameters.fpp:262
$:GPU_DECLARE(create='[p0ref, rho0ref, T0ref, ...]')

fed by bub_pp%p0ref etc. pref/rhoref belong to the Euler-Euler bubble model and the Tait EOS, and the two bubble models are mutually exclusive. Writing pref from the Lagrangian path therefore looks like pref being reached for where p0ref was meant, not a deliberate coupling.

Suggested resolution

Delete the assignment. It has no reachable consumer, the EL model already carries its own reference pressure, and removing it stops a user-supplied case parameter being silently clobbered. If some Lagrangian path genuinely needs a reference pressure, it should use p0ref.

Worth noting separately: an entire equation-of-state model (model_eqns = 4, the Tait/4-equation path) has zero regression coverage. That is the reason this was so hard to reason about — nothing would have caught a change to it either way.

[sbryngelson]

Correction to my previous comment. I listed the hardcoded ICs as the one live read of the global pref. That is wrong — they do not touch the global at all.

Both files declare a local pRef in their variables macro:

#:def Hardcoded2DVariables()
    real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph   ! 2dHardcodedIC.fpp:5

    real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph, Mach   ! 3dHardcodedIC.fpp:3

and assign it before use, in case 204 (Rayleigh-Taylor):

pRef = 1.e5_wp
pInt = pRef
...
q_prim_vf(eqn_idx%E)%sf(i, j, 0) = pref + rhoH*9.81_wp*(1.2_wp - y_cc(j))

Fortran identifiers are case-insensitive, so pRef and pref are the same name. The reference on the last line resolves to the local, which shadows the module global. The hydrostatic initial condition uses a hardcoded 1e5 Pa, exactly as the surrounding hardcoded constants do, and is entirely independent of the case parameter.

Revised conclusion

The global pref has no live readers at all outside the model_eqns = 4 Tait paths, which are never exercised. Same for rhoref. Both become completely dead once the four-equation model is removed, and no hardcoded-IC file needs touching.

That also strengthens the case for deleting m_bubbles_EL.fpp:381 outright: it writes a global that, after the four-equation removal, nothing anywhere reads.

Worth fixing regardless

The spelling flips within a few lines — pRef where it is assigned, pref where it is read. That is legal and correct, but it reads like a reference to the module global, which is precisely how I misread it. Consistent casing would make the shadowing obvious.