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

[qodo-code-review]

ⓘ You've reached your Qodo monthly free-tier limit. Reviews pause until next month — upgrade your plan to continue now, or link your paid account if you already have one.

[sbryngelson]

very impressive! but we're definitely going to need a Riemann solve refactor before merging this. it's all just too cumbersome.

[sbryngelson]

Heads-up on what changed in master since this branch was cut, since the conflicts here are structural rather than mechanical:

1. Unify ICPP STL onto the shared IB model path #1546 unified the flux arrays to x-only: the per-direction flux_rs{x,y,z}_vf arrays are gone; Riemann code now writes flux_rsx_vf(${SF('')}$, ...) with permuted indices. The HLLD hypoelastic anchor block here still emits flux_rs${XYZ}$_vf inside its Fypp direction loop, so the y/z iterations reference arrays that no longer exist — that block needs reworking against the new architecture before a merge with master can compile.
2. Named values for enumerated case parameters and AST-based analytic IC codegen #1550 introduced named parameter values + generated namelists: riemann_solver comparisons in master are now riemann_solver == riemann_solver_hll etc., and the namelist statement in m_start_up comes from a generated include — new parameters get registered in toolchain/mfc/params/definitions.py and appear automatically.

We attempted a maintainer-side conflict resolution and got everything mechanical done, but the HLLD flux block (and a couple of smaller physics choices, e.g. the elastic flux adjustment and the hll_u_interface advection-source semantics) need your judgment. Happy to share the partial resolution as a reference branch if useful — just say the word.

[ChrisZYJ]

Thanks for the update. I'm working on the merge together with a few other changes - will commit soon

[sbryngelson]

Following up on my earlier note about #1546 (x-only flux arrays) — the Riemann refactor I mentioned has now landed (#1550–#1556), and since you said you're mid-merge, here is a map of where your touched regions now live on master so you can relocate rather than re-resolve:

src/simulation/m_riemann_solvers.fpp is now a ~116-line dispatcher. The solver bodies were split out, and your hunks map as follows:

• s_hll_riemann_solver (your interface-consistent HLL changes) → src/simulation/m_riemann_solver_hll.fpp
• s_hllc_riemann_solver (your 12 hunks incl. low-Mach/ADC work) → src/simulation/m_riemann_solver_hllc.fpp
• s_hlld_riemann_solver and your new s_hypo_hlld_riemann_solver (anchored dual-pass) → src/simulation/m_riemann_solver_hlld.fpp
• s_populate_riemann_states_variables_buffers, s_initialize_riemann_solver, s_finalize_riemann_solver, and the viscous source-flux routines → src/simulation/m_riemann_state.fpp; your new s_finalize_nc_iface_vel belongs there next to s_finalize_riemann_solver
• s_riemann_solver" (dispatch), s_initialize/finalize_riemann_solvers_module→ stay inm_riemann_solvers.fpp(wire theriemann_hypo_ADC`/dual-pass dispatch there)

New parameters must go through the registry now. Your hand-edited declarations and m_mpi_proxy broadcasts for riemann_hypo_ADC, ADC_kappa, hll_u_interface, hypo_hll_interface_rhs will conflict: declarations, namelist bindings, and MPI broadcasts are auto-generated from toolchain/mfc/params/definitions.py (_r() + _nv(); physics constraints in case_validator.py). Default assignments stay manual in s_assign_default_values_to_user_inputs; GPU_DECLARE lines for registry scalars live in m_global_parameters_common.fpp (declare directives must be in the declaring module). See docs/documentation/contributing.md. Derived internal flags (hypo_nc_*, adv_src_*, use_nc_iface_vel) can stay as manual declarations in m_global_parameters.fpp.

Other mappings: your m_hypoelastic split (*_finite_diff_per_sweep/*_iface/*_axisym_geom_iface) still targets m_hypoelastic.fpp (unchanged home, but diff against current master — the energy-coupling code there was recently consolidated); m_checker.fpp, m_rhs.fpp, m_time_steppers.fpp, m_variables_conversion.fpp keep their homes. Your SPLIT_LONG_DIRECTIVE macro patch still applies — acc_macros.fpp/omp_macros.fpp retain the $:acc_directive/$:omp_directive emission points.

Two cautions that still stand: (1) the x-only flux architecture concern from my earlier comment — everything writes flux_rsx_vf with permuted indices now, so the per-direction flux_rs${XYZ}$_vf writes in the hypo anchor blocks need rethinking, and that is a design decision I don't want to guess for you; (2) please regenerate your 38 golden files against current master rather than carrying the old ones over — an unexplained golden diff after the refactor would mask exactly the kind of silent regression we're trying to avoid.

Happy to help with the mechanical relocation once you push your merge — ping me here.