Skip to content

Pressure relaxation inverts the stiffened-gas isentrope in closed form; a generic EOS needs a numerical inverse #1785

Description

@sbryngelson

Split out from #1762, which centralized the isentrope parameters but not the inversion. This is
the one place where a second equation of state needs genuinely new numerics rather than a new
coefficient supplier.

What is there now

s_equilibrate_pressure (src/simulation/m_pressure_relaxation.fpp) drives a Newton iteration on
the equilibrium pressure. Each step needs the density on each phase's isentrope at the current
pressure, and the derivative of that density:

rho_K_s(i) = alpha_rho(i)/max(alpha(i), sgm_eps) &
             *((pres_relax + isentrope_B(i))/(pres_K_init(i) + isentrope_B(i)))**(1._wp/isentrope_n(i))
f_pres  = f_pres  + alpha_rho(i)/rho_K_s(i)
df_pres = df_pres - alpha_rho(i)/(isentrope_n(i)*rho_K_s(i)*(pres_relax + isentrope_B(i)))

Both come from the stiffened-gas isentrope p + B = const * rho**n having a closed-form inverse
rho(p) and a closed-form drho/dp. #1762 made isentrope_n and isentrope_B first-class
(f_isentrope_exponent, f_isentrope_pressure, f_pressure_on_isentrope), so the parameters now
have one definition - but the inversion above is still written for that particular functional form.

Why a generic EOS cannot reuse it

For a Mie-Gruneisen or JWL fluid the isentrope has no analytic inverse. JWL's isentrope

$$p_s(\rho) = A e^{-R_1 \rho_0/\rho} + B e^{-R_2 \rho_0/\rho} + C (\rho/\rho_0)^{\omega+1}$$

cannot be solved for rho in closed form, so rho(p) must come from a numerical inversion, and
drho/dp from differentiating the EOS rather than from a power law. That is new solver code, not a
new coefficient table - and it sits inside a Newton loop in a per-cell device routine, so it is a
nested iteration on the GPU with all the convergence and cost questions that implies.

Suggested scope

  • An operator pair f_isentrope_density(pres, pres_0, rho_0, i) and its derivative, with the
    stiffened-gas closed form as the default implementation so nothing changes today.
  • A bracketed Newton (or Ridders) fallback for EOS families with no analytic inverse, with an
    iteration cap and a documented failure mode - it must not silently return a non-converged density
    inside the outer Newton solve.
  • Convergence and cost measured before it is enabled: this is a per-cell device routine in the 6-eq
    path, and a nested iteration there is the kind of change that regresses the hot loop. Pass an eos_state to s_compute_speed_of_sound #1714 was
    abandoned over a 20% regression from a much smaller change to this area.
  • A manufactured test with a known non-analytic isentrope, so the fallback is exercised without
    needing JWL to land first.

Until then, case_validator.py should refuse model_eqns = 3 pressure relaxation for any fluid
whose EOS lacks an analytic isentrope inverse. Blocks the Mie-Gruneisen step of #1638 only for
6-equation cases; the 5-equation Allaire path that #1638 actually targets does not reach this code.

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