Skip to content

Bubble sound speed: alpha -> 1 should assert, not silently produce Inf #1793

Description

@sbryngelson

Summary

The bubble-mixture branch of s_compute_speed_of_sound is an O(alpha) expansion whose alpha is the
dilute bubbly void fraction the user supplies. When alpha approaches one the expansion has left its
regime, and the code responds by silently producing +Inf rather than saying so.

A floor on (1 - alpha) would be the wrong fix - it would let a bad index or an out-of-regime case run
on quietly. What is wanted is a diagnostic that names the actual condition.

The branch is just the liquid sound speed with a void correction

rho under bubbles_euler is the liquid partial density, (1 - alpha) rho_l, and gamma/pi_inf are
the carrier liquid's own coefficients. So

c^2 = K_l / (rho (1 - alpha))  =  K_l / (rho_l (1 - alpha)^2)  =  c_l^2 / (1 - alpha)^2

that is,

c = c_l / (1 - alpha), with c_l = sqrt(K_l/rho_l) the pure carrier-liquid sound speed.

m_qbmm.fpp:775 forms that same c_l directly - f_bulk_modulus(pres, gammas(1), pi_infs(1))*(1 - alf)/rho

  • for Keller-Miksis, and is independent of alpha as a pure-liquid property should be. Verified numerically:
alpha c_qbmm^2 (= c_l^2) c from this branch c_l/(1 - alpha)
0.00 21981.075 148.2602 148.2602
0.04 21981.075 154.4377 154.4377
0.20 21981.075 185.3252 185.3252
0.50 21981.075 296.5203 296.5203

Written this way three things become readable that the current form hides: it is manifestly the O(alpha)
expansion c ~ c_l (1 + alpha); it diverges as alpha -> 1, so a near-one value cannot be a physical
mixture speed; and it is the same quantity m_qbmm computes, differing only by the void correction, so
the two should not be mistaken for competing formulas. (The speed rising with void fraction is correct for
an O(alpha) carrier-phase correction and is why alt_soundspeed exists separately as Wood's law.)

Scope of validity

The expansion is only good in the dilute limit. Outside it neither the Euler-Euler (bubbles_euler) nor
the Euler-Lagrange bubble model is recommended, so a large alpha reaching this line indicates a case
outside the models' intended regime, a wrong index, or both.

What to do

alpha == 1 yields +Inf, which is not < 0, so it passes the mixture_err check and propagates through
sqrt into the wave speeds; alpha > 1 gives a negative c caught only when mixture_err is on (off by
default in the test base config). Either way the failure surfaces far from its cause.

An assertion or diagnostic naming the dilute-limit violation would report the real problem at the point it
occurs. A clamp would not.

Pre-existing on master (src/common/m_variables_conversion.fpp); not introduced by #1762,
which preserves the semantics exactly while centralizing the expression.

Correction to the original report

First filed claiming this division was the mechanism behind NaNs seen on Frontier CCE OpenACC in #1794.
That was wrong. There alpha is nowhere near one: the Lagrange suites run bubbles_lagrange with
bubbles_euler = F and never enter this branch, and the QBMM case has alpha = 4e-2, so the correction is
1/0.96 - which accounts for the measured 1.7% answer shift and nothing more. That NaN mechanism remains
open in #1794.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions