Fix the azimuthal flux divergence and geometric sources in 3D cylindrical coordinates

[aliencaocao]

Fixes #1817.

Changes

All changes are in the θ (z) direction of 3D cylindrical runs (grid_geometry == 3):

1. src/simulation/m_rhs.fpp: the flux divergence is (F(k−½) − F(k+½)) / (r Δθ). The extra factor of the cell-centre velocity u_θ is removed, and the velocity_val private is removed with it.
2. src/simulation/m_rhs.fpp: the four azimuthal volume fraction source terms use 1/(rΔθ) instead of 1/Δθ. The four terms are the per-fluid α trace and the shared velocity trace, each with its Kapila K∇·u part. Their metric matches the flux divergence next to them.
3. src/simulation/m_riemann_solver_hll.fpp, m_riemann_solver_lf.fpp: the θ sweep writes its geometric momentum source instead of the value that the radial sweep leaves in the buffer. Radial momentum gets −(F_θ(ρu_θ) − p̂) = −ρu_θ², with p̂ the solver's own face pressure. Azimuthal momentum gets F_θ(ρu_r) = ρu_r u_θ. All other entries are zero. HLLC writes the same source in the same sweep. HLL and Lax-Friedrichs include pressure in their momentum flux. The new source removes that pressure with the same expression both solvers use in the radial sweep.

The right-hand side applies −(1/r)·½(G(k−½) + G(k+½)) for this source G. The momentum equations therefore receive +ρu_θ²/r (radial) and −ρu_r u_θ/r (azimuthal).

The six-equation branches (θ phasic energy source, HLLC six-equation θ block) have the same omissions. The case validator rejects model_eqns = 3 in 3D cylindrical coordinates. Those branches are left unchanged.

Why this is correct

The conservative form in cylindrical coordinates (vector and tensor divergence in Wikipedia: Del in cylindrical coordinates) is

∂q/∂t + ∂F_x/∂x + (1/r)∂(r F_r)/∂r + (1/r)∂F_θ/∂θ = S
S(ρu_r) = (ρu_θ² + p)/r        S(ρu_θ) = −ρu_r u_θ/r

F_θ already contains u_θ. The azimuthal difference is divided by rΔθ and multiplied by nothing. The non-conservative α∇·u term needs the same 1/(rΔθ). The radial sweep already cancels the p/r part of S(ρu_r). The θ sweep must therefore supply ρu_θ²/r and −ρu_r u_θ/r. HLL and Lax-Friedrichs include the pressure in their momentum flux. Their face pressure (s_M p_R − s_P p_L)/(s_M − s_P) is therefore subtracted. Both solvers use the same expression for the radial sweep.

Before and after

The first two figures show the uniform 40/60 mixture in a swirling pipe flow (HLLC) from the issue. Both cover one cross section of the pipe from t = 0 to 0.1:

• Before: fluid 1 drifts to between 35.55% and 47.87% by t = 0.1.
• After: fluid 1 is 40% in every cell at every saved step, to every written digit.

The figure below shows solid body rotation with HLL, an exact steady state, at t = 1:

• Before: outward flow of up to 0.040 and pressure errors of up to 0.15.
• After: radial velocity at most 1.3 × 10⁻³ and pressure error at most 1.7 × 10⁻⁴, the level of the discretisation error.

Tests and goldens

The existing 3D cylindrical cases vary only along the pipe. Their azimuthal differences vanish, and their goldens stay unchanged. Three cases are added (3D -> Cylindrical -> Azimuthal -> HLL | HLLC | LF). Each case is the standard three-slab pipe with swirl u_θ = 0.1 and an off-axis wedge at a different pressure. The wedge is patch 4, which the Conduction case already uses. The cases share this patch through one dictionary. Pressure and swirl therefore vary with θ. All four patches hold one uniform 50/50 mixture. The cases therefore measure the azimuthal terms alone, whatever a solver does with volume fractions. The goldens were regenerated on CPU with GCC 16.2, in double precision, without case optimization:

• 4C1D77D7 3D -> Cylindrical -> Azimuthal -> HLLC (new)
• 6FEA09D2 3D -> Cylindrical -> Azimuthal -> HLL (new)
• FBD81266 3D -> Cylindrical -> Azimuthal -> LF (new)

#1910 raises the same metric question for the viscous azimuthal flux. That term is computed elsewhere and is not changed here.

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

PR template credit: junegunn

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[github-actions]

Claude Code Review

Head SHA: c8ecef2

Files changed:

• 4
• src/simulation/m_rhs.fpp
• src/simulation/m_riemann_solver_hll.fpp
• src/simulation/m_riemann_solver_lf.fpp
• toolchain/mfc/test/cases.py

Findings:

• medium src/simulation/m_riemann_solver_hll.fpp (also m_riemann_solver_lf.fpp): The new NORM_DIR == 3 block writes flux_gsrc_rsx_vf for cylindrical geometry, but the m_rhs.fpp hunk reads only flux_n(3) and never a geometric-source flux. So the block is only useful if some code outside the diff reads flux_gsrc_rsx_vf for direction 3 (I could not confirm this from the prepared context), and I could not confirm the array is allocated and sized for dir 3. Please confirm the consumer and the allocation. If the consumer is absent, the block is dead code and the azimuthal momentum coupling is still missing. Confidence is moderate.
• medium src/simulation/m_rhs.fpp: The cylindrical z-flux divergence no longer multiplies by q_prim_vf%vf(cont%end + idir) (velocity_val), so the term is now inv_ds*(flux_face1 - flux_face2) with inv_ds = 1/(dz*y_cc). This is a CFD-result change that applies to every equation j = 1..sys_size. The previous velocity factor looks wrong, since flux_n is already a flux. Per CLAUDE.md, PRs that change CFD results need verification. The only coverage added is the new "Azimuthal" golden cases, so please check the new golden files against an independent reference. Also, 1/(dz*y_cc) is only valid if dz is the angular width, and the same y_cc metric change is applied to local_inv_ds in the non-conservative advection terms; please confirm those terms use the same convention. Confidence is moderate.
• medium src/simulation/m_riemann_solver_lf.fpp: The new block uses eqn_idx%mom%beg + 1 and eqn_idx%mom%end as the radial and azimuthal momentum indices. This is correct only when num_dims == 3, which grid_geometry == 3 requires, so please confirm that is enforced in case_validator.py. The block also takes pres_L/pres_R and s_M/s_P from the LF dissipation form (s_P = max(s_L, s_R) + max(c_L, c_R), s_M = -s_P), so (s_M*pres_R - s_P*pres_L)/(s_M - s_P) reduces to (pres_L + pres_R)/2. That is probably intended, but the NORM_DIR == 2 branch has no equivalent evidence in the diff. Confidence is low.
• medium toolchain/mfc/test/cases.py: The new "Azimuthal" stack pushes uniform_mixture and vel(3) = 0.1 over theta_wedge, and the "Conduction" case spreads **theta_wedge into the same dict. In "Conduction", patch 4 keeps vel(3) = 0.0 and pres = 0.5, as before, so it is unchanged. In "Azimuthal", patch 4 pres is overridden to 1.0, and uniform_mixture overrides patch 4's alpha/alpha_rho. The comment says patches 1-3 are untouched. A uniform velocity vel(3) = 0.1 along theta is not a rigid-body rotation, so it is not an equilibrium in cylindrical coordinates, but it is fine for a regression case. Confidence is low.

[codecov]

Codecov Report

❌ Patch coverage is 71.42857% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.96%. Comparing base (47568a0) to head (bc715da).

Files with missing lines	Patch %	Lines
src/simulation/m_rhs.fpp	50.00%	2 Missing ⚠️
src/simulation/m_riemann_solver_hll.fpp	80.00%	0 Missing and 1 partial ⚠️
src/simulation/m_riemann_solver_lf.fpp	80.00%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2051      +/-   ##
==========================================
+ Coverage   61.90%   61.96%   +0.06%     
==========================================
  Files          86       86              
  Lines       22885    22895      +10     
  Branches     3369     3371       +2     
==========================================
+ Hits        14166    14188      +22     
+ Misses       6220     6205      -15     
- Partials     2499     2502       +3     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_riemann_solver_hll.fpp	618	+11
src/simulation/m_riemann_solver_lf.fpp	520	+11
src/simulation/m_rhs.fpp	1967	-1

Directory	Lines	Diff
simulation	28547	+21
total	47606	+21

[wilfonba]

This all appears to be correct. It also seems to fix some of the issues I was seeing with cylindrical coordinates several months ago. Merge away once the CI has run.