HLL axisymmetric: volume fraction solved in conservative form (missing Kapila source) after PR #800

[sbryngelson]

Summary

PR #800 changed one line in the HLL solver's cyl_coord block for the radial (y) direction:

! m_riemann_solvers.fpp, HLL cyl_coord block, NORM_DIR == 2
- flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
+ flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)

This fixed a real issue (HLL was not applying any cylindrical correction to α), but the fix introduced a different error: HLL now solves the volume fraction equation in pure conservative form rather than the quasi-conservative form required by the Kapila/Allaire 5-equation model.

The math

The Kapila volume fraction equation is:

∂α/∂t + ∇·(αu) = α∇·u

In cylindrical coordinates (r, z), the geometric terms cancel:

∇·(αu) - α∇·u  =  u_r ∂α/∂r

So the physically correct equation is the material-derivative form, identical to Cartesian:

∂α/∂t + ∂(αu_r)/∂r = α ∂u_r/∂r

After PR #800, HLL sets flux_gsrc(adv) = flux_rs(adv) and flux_src(adv) = 0, giving:

Term	Value
Flux divergence	-∂(αu_r)/∂r
Geometric correction	-(αu_r)/r
Non-conservative source	0
Net	∂α/∂t + (1/r)∂(rαu_r)/∂r = 0 ✗

This is the purely conservative form with no Kapila source, which will generate spurious pressure oscillations at material interfaces (the original motivation for the quasi-conservative approach).

HLL may appear to work despite this because numerical diffusion masks the interface error, but it is solving the wrong equation.

Correct HLL fix

The right approach mirrors what HLLC already does:

• Keep flux_gsrc(adv) = 0 (geometric terms cancel analytically)
• Supply flux_src(adv) with a velocity estimate so the non-conservative source α ∂u_r/∂r is applied

For HLL, there is no contact wave, but a natural estimate is the HLL-averaged normal velocity:

u_HLL = (s_L * u_R - s_R * u_L) / (s_L - s_R)

Setting flux_src(adv) = u_HLL would make HLL produce:

∂α/∂t = -∂(αu_r)/∂r + α ∂u_r/∂r = -u_r ∂α/∂r   ✓

The cylindrical K·∇·u correction (already present in s_add_directional_advection_source_terms for the alt_soundspeed path) would then also apply correctly for HLL, since it keys off flux_src.

Context

This issue was identified during the analysis for #801. The HLLC solver does not share this problem — its flux_gsrc(adv) = 0 is correct because it compensates via flux_src(adv) = s_S. See the comment at #801 (comment) for the full derivation.

[ChrisZYJ]

@sbryngelson thanks for looking into this! I went through your comment carefully and think #800 is correct and should not be reverted. #800 might not have fixed all the subtle bugs suggested by issue #801, but #800 fixes a critical error in the HLL regularization, and reverting it would make the simulation completely wrong.

To summarize:

• Fix Axisymmetric HLL #794 and Fix multi-comp axisym HLL #800 "fixed" axisym HLL. HLL axisym vs HLLC axisym vs HLLC 3D were visually very different before the fix and visually the same after the fix. This is strong empirical evidence that Fix multi-comp axisym HLL #800 fixed a critical bug correctly.
• It's neither textbook nor straightforward to explain why Fix multi-comp axisym HLL #800 is needed. In short, it's an HLL regularization (compulsory additional numerical diffusion flux) for the non-conservative terms, that must also be applied to the geometric source terms.

I'm working on a (large) HLLD Hypoelasticity PR that should arrive soon, and I'll attach some notes on non-conservative terms that might benefit future developers, as this thing is currently not well-documented and confusing. But for now, here's an AI-generated summary that explains the main ideas:

======

The cancellation argument for the continuous equation is correct — I agree that Dα/Dt = 0 holds in cylindrical coordinates and that HLLC's treatment is sound. However, I believe the conclusion about PR #800 rests on a misidentification of what flux_rs(adv) represents in HLL.

What flux_rs(adv) actually is

MFC's HLL splits the volume fraction numerical flux into two arrays (m_riemann_solvers.fpp, lines 693–699):

flux_rs(adv_i)     = (αL - αR) * s_M * s_P / (s_M - s_P)   ! numerical diffusion
flux_src_rs(adv_i) = (s_M*αR - s_P*αL) / (s_M - s_P)       ! HLL-averaged interface α

flux_rs(adv) is not the full conservative flux αu_r. It is only the numerical diffusion — the Lax-Friedrichs-like dissipation that regularizes α-jumps at material interfaces.

The advective/transport part lives in flux_src_rs(adv), which feeds into a non-conservative source term in m_rhs.fpp (lines 1197–1214):

RHS(α) += v_r * (flux_src_{k-½} - flux_src_{k+½}) / Δr   ≈  -v_r · ∂α/∂r

This already provides quasi-conservative material-derivative transport (Dα/Dt ≈ 0) without any geometric correction — consistent with the cancellation you derived.

flux_src(adv) ≠ 0 for non-MHD HLL

The comment states "HLL: flux_src(adv) = 0". The line that zeros flux_src_rs(adv%beg) (line 762) is inside the if (mhd) block. For standard non-MHD runs, all advection indices carry the nonzero HLL-averaged interface α. The non-conservative source is active.

Why the geometric correction is needed — on the diffusion, not the transport

The cancellation argument applies to the transport operator. Agreed — and that transport is handled by the non-conservative source (term above), which correctly needs no flux_gsrc.

But flux_rs(adv) is numerical diffusion (HLL regularization), not transport. It is the Lax-Friedrichs-type dissipation inherent to HLL — a real numerical flux through cell faces that smooths α-discontinuities. In cylindrical coordinates, any face-flux F requires the full cylindrical divergence:

(1/r) ∂(rF)/∂r  =  ∂F/∂r  +  F/r
                    ─────      ───
                    term (1)   term (2) = flux_gsrc

Before PR #800, only term (1) was applied to this diffusion flux — Cartesian weighting regardless of radius. PR #800 adds term (2), completing the cylindrical divergence for the diffusion operator.

Why HLLC doesn't need this

(Chris: HLLC uses a different non-conservative regularizing flux from HLL. Again, I'll explain it in my notes)

Why PR #794 worked for single-component but not multi-component

• Single-component: α = 1 everywhere → flux_rs(adv) = (1-1)*... = 0 → geometric correction is zero regardless.
• Multi-component: α jumps at interfaces → flux_rs(adv) ≠ 0 → Cartesian-only diffusion is geometrically inconsistent → instability.

Summary

Term	What it does	Needs geometry?
flux_rs(adv) = HLL diffusion	Regularizes α-jumps	Yes — it's a face-flux
flux_src_rs(adv) = interface α	Quasi-conservative transport via v_r·∂S/∂r	No — cancellation applies

PR #800 applies cylindrical geometry to the diffusion operator, not to the transport. The cancellation argument is correct for transport but does not apply to diffusion.