Give HLL's viscous energy flux a face velocity

[aliencaocao]

Fixes #2046.

Change

In src/simulation/m_riemann_solver_hll.fpp, the loop that exports the HLL interface velocity for the non-conservative terms also writes it to vel_src_rsx_vf when viscous = T. That array is the one the shared viscous routine reads for the τ·u work term. The velocity is the existing HLL trace. It is upwind outside the wave fan and averaged inside it:

Ψ_u = u_L                                   if 0 ≤ s_L
Ψ_u = (s_R u_L − s_L u_R) / (s_R − s_L)     if s_L ≤ 0 ≤ s_R
Ψ_u = u_R                                   if s_R ≤ 0

The formula holds for every velocity component. The condition that selected this loop is kept in one logical, need_nc_vel. The non-conservative export is therefore unchanged.

Why this is correct

MFC's energy equation carries the viscous work in its flux, ∇·[(ρE + p)u − τ·u] (equations.md, Section 2.1). With a consistent face velocity, the internal energy then gains exactly the dissipation Φ = τ:∇u ≥ 0 (Mei, MIT 1.63J lecture notes, Eqs. 4.1.16 and 4.1.18). Ψ_u is the HLL flux of the quantity whose "flux" is u. This is the interface velocity that MFC already uses for HLL's volume fraction and Kapila terms (HLL_HLLC_non_conservative_terms_derivations.md, Section 3). HLLC uses its own star velocity for the same work term. For symmetric wave speeds, it reduces to the arithmetic mean of u_L and u_R. That mean is the face velocity of a central viscous flux (Coralic and Colonius 2014, Section 3.2).

Before and after

The figure below shows the Gaussian vortex from the issue on 64², 128² and 256² cells. Each panel shows the heating rate over the first 0.005 time units against the analytic dissipation Φ = 4μA²(r/R)⁴ exp(−2r²/R²).

• Before: the core overheats, a ring at r ≈ 0.4 to 0.6 cools, and the relative RMS error is 112.9% to 113.9% on every grid.
• After: the heating ring of Φ appears, and no cell cools on any grid (minimum rate approximately −3 × 10⁻⁷, round off). The relative RMS error is 3.8%, 1.7%, 1.5% on 64², 128², 256².

Tests and goldens

The suite's Viscous cases all run the default HLLC solver. The one suite case that combines HLL with viscosity is the viscous_shock_tube example. The golden was regenerated on CPU with GCC 16.2 in double precision without case optimization:

• 70EC99CE 2D -> Example -> viscous_shock_tube

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: 7abeb33

Files changed:

• 1
• src/simulation/m_riemann_solver_hll.fpp

Findings:
No issues found. Checked in step 4: src/simulation/m_riemann_solver_hll.fpp (against m_riemann_state.fpp consumers of vel_src_rsx_vf and the HLLC writers).

[codecov]

Codecov Report

❌ Patch coverage is 28.57143% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.89%. Comparing base (47568a0) to head (7abeb33).

Files with missing lines	Patch %	Lines
src/simulation/m_riemann_solver_hll.fpp	28.57%	3 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2054      +/-   ##
==========================================
- Coverage   61.90%   61.89%   -0.01%     
==========================================
  Files          86       86              
  Lines       22885    22888       +3     
  Branches     3369     3371       +2     
==========================================
+ Hits        14166    14167       +1     
- Misses       6220     6221       +1     
- Partials     2499     2500       +1     

☔ 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	610	+3

Directory	Lines	Diff
simulation	28529	+3
total	47588	+3

[wilfonba]

I ran and visualized the 2D viscous shock tube, and this all looks to be correct. Merge away

[sbryngelson]

can you better explain the code changes and side effects please?

[aliencaocao]

1. The face velocity is now also written for viscous runs. HLL already computed an interface velocity Ψ_u at each face: upwind outside the wave fan and the HLL average inside it. It wrote it only to nc_iface_vel_rsx_vf, and only when the non-conservative terms need it (hypo_nc_mode_interface, or alt_soundspeed without hll_u_interface). The same loop now runs when viscous = T too, and it writes Ψ_u to vel_src_rsx_vf. That is the array s_compute_viscous_source_flux reads for the −τ·u work term in the energy flux, in both the Cartesian and the cylindrical routines. Before this PR, HLL never wrote vel_src_rsx_vf, so the work term multiplied τ by garbage values. HLLC fills the same array with its star velocity, so HLLC was never affected.

2. Refactor so the formula is written once. Ψ_u is computed once into a private scalar vel_interface. That scalar is then stored to one or both arrays. The old condition moves to a logical, need_nc_vel, which is set once before the loop and passed in as firstprivate. Because of this, the non-conservative export uses the same formula on the same inputs, so its values are unchanged.

Side effects (not really sure what you mean by this, so I just listed what changes it makes to the results and runtime)

• Corrects the numerics of HLL + Viscous = T cases, only 1 now in the test suite.
• No memory usage change. vel_src_rsx_vf is already allocated for every solver in m_riemann_solvers.fpp.
• Cost: for viscous HLL runs, one extra store per face per dimension. This is necessary for correctness.