Repository navigation
Lagrange bubble initialization silently overwrites the user-supplied global pref #1710
Description
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 ifThere 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
-
prefis a user input. It is registered attoolchain/mfc/params/definitions.py:739, described as "Reference pressure", documented incase.md(line 137, and line 876 as "Reference pressure for bubble models"), and required bymodel_eqns = 4(definitions.py:405). Overwriting it discards a value the user set. -
It is per-bubble, in a loop.
s_add_bubblesis 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_pressureuses 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_UPDATElist atm_global_parameters.fpp:916, so the overwritten value reaches GPU kernels. -
It is conditional on surface tension. The assignment only happens when
1/Webis not approximately zero, so whether a user'sprefsurvives 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.
Activity
sbryngelson commentedon Aug 8, 2026
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
m_global_parameters_common.fpp:378dflt_realpre_process/m_global_parameters.fpp:520bubbles_euler .and. polytropic->1._wppost_process/m_global_parameters.fpp:419simulation/m_bubbles_EL.fpp:381bubbles_lagrange+ surface tension -> last bubble'sgas_pReads
m_variables_conversion.fpp:105model_eqns == 4(Tait inversion)pre_process/m_assign_variables.fpp:409elseofif (model_eqns /= model_eqns_4eq)— somodel_eqns == 4pre_process/m_icpp_patches.fpp:646model_eqns == model_eqns_4eq2dHardcodedIC.fpp:115,121,3dHardcodedIC.fpp:158,164Consequence
Every EOS-side consumer of
prefis gated onmodel_eqns = 4, andmodel_eqns = 4is never exercised: the test suite only ever setsmodel_eqnsto 2 or 3, and no example or benchmark sets 4 either. No case anywhere setspref.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:381a 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 stillGPU_UPDATEd to device.The likely explanation
The Euler-Lagrange model has its own reference set, declared separately from
pref/rhoref:fed by
bub_pp%p0refetc.pref/rhorefbelong to the Euler-Euler bubble model and the Tait EOS, and the two bubble models are mutually exclusive. Writingpreffrom the Lagrangian path therefore looks likeprefbeing reached for wherep0refwas 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 commentedon Aug 8, 2026
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
pRefin their variables macro:#:def Hardcoded2DVariables() real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph ! 2dHardcodedIC.fpp:5and assign it before use, in case 204 (Rayleigh-Taylor):
Fortran identifiers are case-insensitive, so
pRefandprefare 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
prefhas no live readers at all outside themodel_eqns = 4Tait paths, which are never exercised. Same forrhoref. 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:381outright: it writes a global that, after the four-equation removal, nothing anywhere reads.Worth fixing regardless
The spelling flips within a few lines —
pRefwhere it is assigned,prefwhere 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.