CCE OpenACC faults when the acoustic loop calls s_compute_speed_of_sound

[sbryngelson]

Summary

m_acoustic_src computes its own sound speed instead of calling s_compute_speed_of_sound. It is the
only site in the solver that does so. Under bubbles_euler the two disagree, which means the acoustic
source injects a wave at a different speed than the wave speeds the Riemann solver computes on the same
cells.

Attempting to unify them destabilizes the solver on Frontier CCE OpenACC, so the disagreement needs a
physics ruling before either form can be called correct. This issue records the derivation and the
measurements so none of it has to be rediscovered.

The two expressions are the same EOS

m_acoustic_src (src/simulation/m_acoustic_src.fpp) forms the mixture coefficients and then writes the
Tait form by hand:

call s_compute_mixture_coefficients(myalpha_rho, myalpha, myRho, small_gamma, B_tait, qv_dummy)
small_gamma = 1._wp/small_gamma + 1._wp
c = sqrt(small_gamma*(pres + ((small_gamma - 1._wp)/small_gamma)*B_tait)/myRho)

With n = 1/Gamma + 1, (n - 1)/n = 1/(1 + Gamma), so this is exactly
sqrt(((Gamma + 1)p + Pi)/(Gamma rho)) - i.e. f_bulk_modulus(p, Gamma, Pi)/rho, the plain branch of
s_compute_speed_of_sound. The downstream c/(small_gamma - 1) is likewise a pure rename: n - 1 = 1/Gamma,
so it equals c*Gamma.

So there is no EOS disagreement. The entire behavioural delta is what s_compute_speed_of_sound adds
on top of the plain branch and the hand-rolled code never applied:

• the bubble dilution c = c/(1 - adv(num_fluids)) under model_eqns_5eq .and. bubbles_euler
• the alt_soundspeed Wood's-law harmonic mean
• the model_eqns_6eq arithmetic mean

The open question is therefore specifically whether the dilution closure belongs in the acoustic source
term, not which EOS is right.

Evidence

Switching m_acoustic_src to s_compute_speed_of_sound:

• moved the bubbly QBMM golden CE9DBA3F by 1.7% on AMD MI210 / amdflang OpenMP offload, which is the
  1/(1 - alpha) dilution to within rounding (alpha = 4e-2 predicts ~2.1%). Non-bubbly acoustic goldens
  were unaffected.
• NaN'd at step 50 on Frontier CCE OpenACC in all three suites that combine an acoustic source with a
  bubbly medium: B9553426 and 4A1BD9B8 (2D Lagrange bubbles, acoustic_source = T, mag = 2e04) and
  CE9DBA3F (1D QBMM bubbles). Reverting the change returns all three to green - 41 tests pass, 0 fail,
  including the full Lagrange suite and a 28-case CCE-critical sweep.

The NaN mechanism is still open

An earlier version of this report attributed the NaN to the 1/(1 - alpha) dilution. That was wrong.
In MFC a 1/(1 - alpha) always denotes the dilute bubbly void fraction, and in these cases it is small:
the two Lagrange suites run bubbles_lagrange with bubbles_euler = F and so never enter that branch,
and the QBMM case has alpha = 4e-2, giving 1/0.96. That accounts for the 1.7% answer shift and nothing
more - it cannot produce Inf or NaN.

With model_eqns = 2, alt_soundspeed = F and mixture_err = F (the test base config), the new call
reduces for the Lagrange cases to the identical plain branch, f_bulk_modulus(p, Gamma, Pi)/rho. So the
arithmetic delta there is only floating-point re-association, which should not destabilize anything.

The remaining common factor across all three failures is structural rather than numerical: the change
passes myalpha - an array whose bound num_fluids is a device global - into s_compute_speed_of_sound,
a GPU_ROUTINE(parallelism='[seq]'), from inside the acoustic parallel loop. That is the same construct
class that has twice caused CCE-only failures in this codebase, once as a compile error
(ftn-7066 ... Global in accelerator routine without declare) and once as a runtime GPU memory fault
with a clean build (the m_thermochem case). Clean build plus runtime NaN on CCE only, with amdflang
OpenMP offload green on the same tree, fits that pattern.

This is a hypothesis, not a measurement. Discriminating it needs a variant that computes the mixture sound
speed without passing an array into a device routine from that loop.

What a fix needs

1. A ruling on whether an acoustic source in a bubbly medium should inject at the diluted mixture speed or
   at the carrier-liquid speed. Note m_qbmm.fpp:775 deliberately goes the other way - it multiplies by
   (1 - alf) to recover the pure liquid speed from the partial density, which is correct for
   Keller-Miksis. These are different physical quantities and neither should be folded into the other.
2. If the dilution does belong, Bubble sound speed: alpha -> 1 should assert, not silently produce Inf #1793 has to be fixed first or the source term is unstable wherever the
   void fraction approaches one.
3. Frontier CCE OpenACC validation of the bubbly acoustic cases before merge, not after.

Reverted in #1762 as e07ee640; the rest of that PR stands and is CCE-validated without it.

[sbryngelson]

Result: a variant that works, mechanism still NOT isolated

Frontier CCE OpenACC, job 5381470, at branch head plus the patch below: clean compile (372 s, zero ftn-
diagnostics), 3 passed / 0 failed on 4A1BD9B8, B9553426, CE9DBA3F.

The variant applies the reverted acoustic change but computes the sound speed inline from scalars:

call s_compute_mixture_coefficients(myalpha_rho, myalpha, myRho, gamma_mix, pi_inf_mix, qv_dummy)
c = sqrt(f_bulk_modulus(q_prim_vf(eqn_idx%E)%sf(j, k, l), gamma_mix, pi_inf_mix)/myRho)

instead of call s_compute_speed_of_sound(..., myalpha, c). Note this computes only the plain branch, so
it omits the bubble dilution - CE9DBA3F is a control here, tracking the restored golden, not a test of
the dilution closure.

Correction: "passing an array into a [seq] routine" is NOT the mechanism

An earlier comment here proposed that the failure came from passing myalpha - an array bound by the
device global num_fluids - into a GPU_ROUTINE(parallelism='[seq]') from that loop body. That is
wrong, and there is a live counterexample one line above the change. s_compute_mixture_coefficients
takes alpha_rho_K, alpha_K as dimension(num_fluids) and is called from the same loop body in the
reverted baseline that works and in the patched version that works. So the array crossing cannot be the
mechanism on its own.

Two follow-on candidates were checked and are also eliminated by counterexample:

candidate	eliminated by
array bound by a device global crossing into a [seq] routine	s_compute_mixture_coefficients does exactly this in both working versions
an inner GPU_LOOP(parallelism='[seq]') over num_fluids inside the callee	s_compute_mixture_coefficients contains one as well
cray_inline=True / function_name= present on one decorator and not the other	ACC_ROUTINE (src/common/include/acc_macros.fpp) ignores function_name and builds '!$acc routine ' + clause_val; cray_inline reaches !DIR$ INLINE only on Cray CPU builds, so under MFC_OpenACC both routines emit an identical !$acc routine seq

What is actually established

Only this: that specific call to s_compute_speed_of_sound from the acoustic loop had to go. The
mechanism is not identified. Do not generalize this to "pass scalars, not arrays" - that would send
someone refactoring call sites that were never implicated, including ones this very experiment proves are
fine.

Remaining untested candidates, in rough order of plausibility:

1. A second [seq] device-routine call in the same loop body (call count or nesting, not the argument
   shape).
2. Something in s_compute_speed_of_sound's own body - it reads considerably more device globals than
   s_compute_mixture_coefficients (chemistry, relativity, alt_soundspeed, model_eqns, mpp_lim,
   mixture_err, sgm_eps, gammas, pi_infs) across deeper branching.

Isolating it needs one more variant - e.g. calling a stripped copy of s_compute_speed_of_sound with only
the plain branch, which would separate (1) from (2). Not run: the branch has a green path without it and
the allocation is better spent elsewhere. Anyone picking this up has the experiment pre-designed.

[sbryngelson]

Physics resolved; this issue now tracks the CCE mechanism only

The dilution question is settled and shipped in #1762 as 92091a23: an acoustic source in a
bubbly medium injects at the mixture speed. Since alf is dilute by construction the diluted and
undiluted forms are close, so adopting the shared rule costs little and removes the disagreement with the
wave speeds the Riemann solver computes on the same cells. 39 bubbles_euler goldens move by ~1.7% at the
source cell.

Verified at 92091a23 across every case carrying an acoustic source or bubbles - 119 requested, 119
selected, 119 passed - on both backends: Frontier CCE OpenACC (job on frontier08683, clean compile,
zero ftn- diagnostics) and MI210 amdflang OpenMP offload. That is also the first CCE run to exercise the
dilution branch at all; every earlier green run either omitted it or reverted it.

What remains

The correction is inlined behind the same guard rather than calling s_compute_speed_of_sound,
because that call still faults CCE OpenACC from this loop. Inlining demonstrably avoids whatever the call
was hitting, but the mechanism is not identified, so the duplication of a three-line guard stands as a
deliberate TODO in the source comment. Fold it back into the call once the mechanism is known.

Re-stating the constraint on any future write-up, since it is easy to get wrong: this is not "arrays
must not cross into device routines." s_compute_mixture_coefficients passes myalpha_rho, myalpha -
dimension(num_fluids), a device-global bound - into a [seq] routine one line above, in every version
that works, including the shipped one. The eliminated theories and the untested candidates are listed in
the earlier comment.