Anchored Dual-pass HLLD for Hypoelasticity (+ HLLC and interface-consistent HLL)

[ChrisZYJ]

Description

Adds:

1. Hypoelasticity: Anchored Dual-pass HLLD
2. Hypoelasticity: HLLC
3. Hypoelasticity: HLL option (interface-consistent)
4. HLL option (alpha div U) so non-conservative treatment aligns with HLLC

Key Design Choices

Separate HLLD Riemann Solvers

At a glance it might be tempting to combine HLLD MHD with dual-pass hypoelasticity HLLD, but keeping them separate makes the code cleaner and much easier to maintain because:

1. Unlike HLL or HLLC, HLLD is a class of HLLD-type solvers, with formulas and states dependent on the eigenstructure of the governing equations, so the inner states' equations are completely different for MHD vs Hypoelasticity.
2. HLLD hypoelasticity has a newly developed dual-pass anchored form, making it different from any convenional HLLD Riemann solver. The anchored forms are necessary for the non-conservative hypoelasticity terms, which MHD does not have.
3. MHD and Hypoelasticity deal with completely different physical regimes with different governing equations, and any changes or new physical models added in the future will not apply to both modules at once.

Riemann Source Terms

For the non-conservative terms, unlike the usual governing equations that only need div U i.e. du/dx, dv/dy, dw/dz (alpha div U, K div U, etc.), Hypoelasticity has cross terms like du/dy, so we must also pass those Riemann-consistent traces from Riemann solver to the rhs. (The old Hypoelasticity code with the HLL Riemann solver uses finite difference for non-conservative rhs, which provides enough stability given that HLL smears the interface immediately, so there wasn't a need to pass the du/dy traces before this PR. But that does not work for HLLC/HLLD for Hypoelasticity.)

Also grouped/named the condition branches (with lots of comments within the code):

Branch	Face quantity read	RHS formula per $\alpha_k$	K*div(u) velocity source
adv_src_alpha_iface	flux_src_n(dir)%vf(j_adv) = per-fluid $\Psi_{\alpha_k}$	$u_\text{cell} \cdot \Delta\Psi_\alpha / \Delta x$	nc_iface_vel_n(dir)%vf(1)
adv_src_vel_iface	flux_src_n(dir)%vf(adv\%beg) = shared $\Psi_u$	$\alpha_k \cdot \Delta\Psi_u / \Delta x$	Same flux_src_n slot (already $\Psi_u$)
adv_src_none	—	Skipped (HLLD handles internally)	—

The derivations, meanings, and usage of the Riemann source variables are not straightforward. I've added some hopefully very helpful notes in misc/dev_notes for future developers (or AI agents; directing them to my notes should help them make fewer mistakes with the source terms) in terms of the understanding and derivations for the HLL/HLLC non-conservative fluxes, and their variable mapping for Riemann solvers and RHS.

Backwards Compatibiilty

• All default behaviors preserved exactly (newly added features as options)
  • The only exception is the removal of an incorrect ad-hoc fluids-limit guard that affects only Hypoelasticity HLL
• All existing usage of Riemann and rhs source terms are preserved. No refactor is done to keep the scope of this PR limited (any refactoring would touch most of the existing HLLC functionalities)

Type of change

• New feature

Testing

• All tests passed locally on CPU and Nvidia GPU, and on Frontier.

• Smooth Eigenmode Convergence

• Weak Solution Comparison (Rodriguez & Johnsen (2019) §5.3(b))

• Weak Scaling on Frontier

Checklist

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed
• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

AI code reviews

Reviews are not triggered automatically. To request a review, comment on the PR:

• @coderabbitai review — incremental review (new changes only)
• @coderabbitai full review — full review from scratch
• /review — Qodo review
• /improve — Qodo code suggestions
• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Add label claude-full-review — Claude full review via label

[sbryngelson]

One note on CI. The last two runs each failed a single Frontier job at the Fetch Dependencies step, before anything gets built or tested (gpu-omp 1/2, then gpu-acc 2/2). Both failures show the same error on the same runner, frontier-5:

error: Failed to install: pandas-3.0.3-cp313-cp313-manylinux_2_24_x86_64.manylinux_2_28_x86_64.whl
  Caused by: failed to open file ``/tmp/uv-cache-sbryngelson/archive-v0/h07_8qgFw9kMCAsSlH4jZ/pandas-3.0.3.dist-info/METADATA``: No such file or directory

Maybe a corrupted uv cache entry on that runner's node? Could you help me re-run the one failed job? All the tests should pass.

#1630

[sbryngelson]

Hi @ChrisZYJ — #1572 (Phase-3 Riemann hot-path decomposition) restructured/removed the exact solver code this PR builds on, which is what put m_riemann_solver_hll/hllc/lf.fpp and m_riemann_state.fpp into conflict.

I resolved the merge of current master into this branch for you. It's pushed to:

sbryngelson:pr1414-hypo-hlld-master-merge

Its first parent is your current head (815f1e2f4), so you can take it as a clean fast-forward if you want — no obligation:

git fetch https://github.com/sbryngelson/MFC.git pr1414-hypo-hlld-master-merge
git checkout hypo_hlld
git merge --ff-only FETCH_HEAD
git push

(Or just git pull https://github.com/sbryngelson/MFC.git pr1414-hypo-hlld-master-merge.) It's your PR, so please review the resolution before adopting it — the choices below involve your physics, and I have not compiled or run it (fypp preprocesses cleanly and the generated Fortran is block-balanced, but goldens are untested).

The core tension

#1572 deleted the HLLC/LF hypoelastic-energy blocks as dead code, justified by case_validator.py restricting hypoelasticity to HLL (riemann_solver == 1). This PR lifts that restriction to [1, 2, 4] (HLL/HLLC/HLLD), so those blocks are live again. The resolution keeps your feature and layers master's refactors underneath.

Choices made (10 conflict hunks)

File / region	Resolution	Rationale
LF hypoelastic block	took master (deleted)	LF isn't in [1,2,4], so the block is dead on both sides; merged LF is identical to master
HLL inline energy macro	took master's helper	s_compute_hypoelastic_interface_energy (called just after the loop) already reproduces your inline physics including the shear_indices doubling — keeping the inline macro too would double-count
HLLC declarations	took yours (superset)	your ADC/star-state locals are needed downstream
HLLC 6-eq & 5-eq private lists	took yours, and added master's alpha_lim_L/R to the 5-eq list	GPU-safe superset; the merged 5-eq body uses both your ADC vars and master's alpha_lim arrays
HLLC 6-eq & 5-eq hypoelastic energy	took yours (live feature)	validator now permits HLLC hypoelasticity
HLLC 5-eq mixture loads	took master (alpha_lim_L/R = ...)	the auto-merged helper call accumulates mixture props from alpha_lim; your inline accumulation would double-count against it
HLLC flux else (not-hypoelasticity)	took yours	your else holds the full base momentum+energy+elastic flux; master's conflict-side tail alone would leave momentum/energy flux uncomputed
m_riemann_state	union	kept master's three shared helpers and your templated s_finalize_riemann_solver (+_hatR); dropped master's now-redundant plain s_finalize_riemann_solver decl since your template supersedes it

Two things worth your eyes

1. hypo_energy_guard gating. Your HLLC path uses @:compute_hypo_elastic_energy (which honors hypo_energy_guard), while HLL now goes through master's helper (which hard-gates on G > verysmall). Both are live, and I preserved each side's own choice — but HLL and HLLC now treat the soft-G case differently. You may want them consistent.
2. Cosmetic: an unused Re_max local survives in HLLC (Phase-3: Riemann hot-path decomposition into shared GPU device helpers #1572 had removed it). Harmless; drop at your leisure.

Happy to adjust any of these if you'd prefer a different call. cc @sbryngelson

[ChrisZYJ]

@sbryngelson Thanks again for preparing the master merge. I reviewed and adopted the resolution, retained the requested refactor and cleanup, and addressed the remaining energy and validation consistency items, including the permanent per-side near-zero-G guard and thermal/kinetic enthalpy for the base acoustic sound speed.

All CI tests have passed. Could you please take a final look and approve/merge if satisfied? Thanks!

[sbryngelson]

bumping this so I get to it sooner than later.

[sbryngelson]

Three small robustness items in the anchored HLLD path, from working through the
solver algebra. None of them affects the results in this PR as it stands — the
first two are latent and the third is a guard against future regressions.

1. The inner wave speeds are not clipped into the outer fan.

In m_riemann_solver_hypo_hlld.fpp the inner speeds are set as

sqrtC_NC = sqrt(C_NC)
S_Lstar = S_M - sqrtC_NC/rhoL_star
S_Rstar = S_M + sqrtC_NC/rhoR_star

and the only degeneracy test beforehand is

degenerate = (abs(denomA) < verysmall .or. &
              .not. (S_L - verysmall <= S_M .and. S_M <= S_R + verysmall))

which checks the contact against the outer fan but never the inner waves against
it. C_NC is built from the anchor impedance while rhoL_star/rhoR_star come
from the reconstructed face state, so where the anchor cell and the face state
hold different materials — i.e. at every one-cell interface — the ratio can
exceed c_f. Sweeping the material table, a water face state anchored on a
stiffer solid gives sqrtC_NC/rho_star ≈ 2670 m/s against c_f ≈ 1625 m/s, so
S_Lstar < S_L.

The failure is quiet rather than loud. f_hlld_wave_zone classifies the
interface with an ordered if-chain (0 <= S_L, then 0 <= S_Lstar, …), so an
inverted ordering just makes the outer-star wedge unreachable and the flux is
taken from the next wedge. No NaN, no crash, just a flux from the wrong region of
the fan. Clamping S_Lstar below by S_L and S_Rstar above by S_R fixes it
in two lines, and it subsumes the C_NC -> 0 repair, since collapsing onto S_M
is the limiting case of the clamp.

None of the material pairs in the current example cases trigger this, so it is
latent, not active.

2. The degenerate-shear guard is dimensional.

C_NC = rho_hat*(G_hat + tau_nn_hat)
if (C_NC < verysmall) then

C_NC has units of kg·Pa·m⁻³ and is ~1.8e12 for epoxy, so comparing it against
an absolute floating-point tolerance carries no fixed relation to the material.
The risk is that the guard is too permissive rather than too strict: a
near-fluid cell from volume-fraction undershoot can sit at C_NC ~ 1e-6 in SI,
far above verysmall, so the guard stays off while 1/sqrtC_NC ~ 1e3 amplifies
the double-star states, which go as dtau/(2*sqrtC_NC). Comparing against the
acoustic modulus instead, e.g. C_NC < eps*rho_hat*gamma_hat*(p_hat + pi_hat),
or gating on c_s/c_f, makes it scale-free. One line.

3. A regression test with a sharp criterion.

The reason HLLD beats HLLC at a shear-bearing contact is that its momentum and
energy fluxes carry the same tangential star state — F[rho v] = -tau** and
F[rho E] = -v** tau** — so the kinetic term cancels the velocity-linear part of
the energy update identically. HLLC's carry different weights, (1-beta) in the
momentum flux against beta in the energy flux, and the mismatch is
first order in the shear jump.

That gives a test with a much sharper pass criterion than an eyeball comparison:
set up a static contact (u = 0, tau_xx = 0, p_T uniform, tau_xy jumping),
where the exact solution has u = 0 and p = p_0 so any pressure is spurious,
then sweep the shear jump over a few doublings and fit the order. HLLD should
come out at 2, HLLC at 1. In a 1D reference implementation the fitted orders are
1.99–2.00 and 1.00–1.02 respectively, clean enough that asserting 2.0 ± 0.1 for
HLLD would catch any future break of that flux pairing — a refactor, a GPU port,
a new reconstruction — that a smoke test would miss.

Related, for the docs: since the advantage is this pairing, riemann_hypo_ADC
blending toward HLL trades it away. Probably worth a sentence at the option so
that is explicit rather than implicit.

[ChrisZYJ]

I've made the changes as requested:

1. Inner-fan ordering

The raw S_Lstar and S_Rstar values are now checked against [S_L, S_R] for each anchor. I initially tried clamping them, but a wave speed and its star-state jump are coupled through the Rankine–Hugoniot relation. Changing only the speed would make that wave internally inconsistent.

If the five-wave ordering is invalid, the complete anchored contribution now uses the existing HLL fallback, or the corresponding one-sided physical flux when the outer fan does not straddle zero. The exported face state and NC interface velocity use the same fallback-safe wave zone, so they cannot sample the invalid five-wave fan. The valid HLLD path is unchanged, while the exact or cancellation-level zero-shear case remains a separate finite HLLC-like collapse.

Added regression test C531DC93 - creates an anchor/face mismatch that distinguishes the two policies. HLL fallback gives a measured one-step transverse-momentum $L_\infty$ change of 0.001; evaluating the rejected speed-clipped five-wave fold gives approximately 0.0084515.

2. Dimensionless shear-degeneracy guard

I replaced the former absolute C_NC threshold with a cancellation-relative criterion:

G_eff   = G_hat + tau_nn_hat
delta_G = 4*epsilon(1._wp)*(abs(G_hat) + abs(tau_nn_hat))

When the scale is nonzero, this is equivalent to testing the dimensionless signed cancellation ratio G_eff/(abs(G_hat) + abs(tau_nn_hat)) against ±4*epsilon(1._wp). The result is independent of the units used for stress.

The solver classifies the state as follows:

• Resolved positive G_eff, with positive representable rho_hat*G_eff: retain the HLLD shear subfan, subject to the independent fan-ordering check.
• -delta_G <= G_eff <= delta_G: use the finite HLLC-like zero-shear limit.
• Materially negative or NaN G_eff: reject the non-real shear subfan and use HLL fallback.
• Positive G_eff for which rho_hat*G_eff underflows to zero: use the numerical zero-shear limit; nonpositive or NaN density falls back.

The factor four is a small conservative allowance for the two mixture products, their accumulation into G_hat, and the final signed sum with tau_nn_hat.

I did not use an acoustic-relative cutoff. This guard is intended to detect cancellation and loss of floating-point representability, not to classify the material. Scaling by the acoustic modulus would eventually classify every sufficiently soft but physically valid solid as fluid-like and suppress a resolved shear subfan, even though a small positive shear stiffness is not singular.

Volume-fraction admissibility also remains separate. MFC already handles reconstructed-state admissibility through mpp_lim; solver-local clipping or an additional ad hoc alpha threshold would modify diffuse-interface states and could compromise consistency and accuracy.

Added 2 regression tests to cover the changed outcomes:

• 10CCE589 verifies that a small but resolved positive stiffness retains HLLD.
• A96DA754 verifies HLL fallback for materially negative effective stiffness.

3. Shear-contact regression

I added Convergence -> HypoShearContact -> {HLLD, HLLC}. It holds the grid fixed, sets the initial normal velocity to zero, retains a nonzero uniform tangential velocity, and sweeps the tau_xy jump over five amplitudes spanning four doublings. The diagnostic measures the baseline-subtracted pressure response after one step. The exact tangential evolution preserves $u=0$ and $p=p_0$, so this pressure response is purely numerical.

The fitted orders are:

• HLLD: 2.000 (required 2.0 ± 0.1)
• HLLC: 0.998 (required 1.0 ± 0.2)

The HLLC leg is intentionally asserted as a control. If the uniform tangential velocity were accidentally removed, both solvers would appear quadratic and the HLLD-only test would lose its discriminating power.

4. ADC documentation

The documentation now states explicitly that ADC trades contact and shear resolution for HLL robustness and is default-off/recommended-off.

I also corrected the documented ADC_kappa direction. It scales the reference jump in the sensor denominator, so smaller values produce more HLL blending, while larger values keep the flux closer to HLLC/HLLD.

[ChrisZYJ]

The one failed documentations CI is unrelated to this PR - it's Lychee link checking for a file unchanged during this PR