Skip to content

[WIP] Fix the HLLC scheme on moving grids and its Jacobians - #2943

Open
rois1995 wants to merge 6 commits into
fix_underrelaxation_factorfrom
fix_hllc_moving_grids
Open

rois1995 wants to merge 6 commits into
fix_underrelaxation_factorfrom
fix_hllc_moving_grids

Conversation

@rois1995

@rois1995 rois1995 commented Sep 30, 2026 •

Copy link
Copy Markdown
Contributor

Proposed Changes

Moving grids. HLLC must transform consistently when the mesh and flow receive the same uniform velocity boost. Some ALE terms were missing or inconsistent, producing frame-dependent fluxes. This PR evaluates the scheme in the moving-face frame and converts the flux back. The comparisons and separate wall/outlet findings are documented in the comments.

The cause is that some terms of the moving face were missing or inconsistent. The sound speeds were shifted by the grid speed, the supersonic branches missed the pressure work of the moving face, and the star-state energy missed the term of the moving face:

SoundSpeed_i -= ProjInterfaceVel;
SoundSpeed_j += ProjInterfaceVel;

Flux[nVar-1] = Enthalpy_i * Density_i * ProjVelocity_i;

IntermediateState[nVar-1] = rhoSL * Density_i * Energy_i - (Pressure_i * ProjVelocity_i - pStar * sM) / (sL - sM);

The flux is now computed as seen from the moving face, where the scheme is the same as on a fixed grid, and then converted back.

Jacobians (the implicit part).

  • For general (non-ideal) gases, one derivative had the wrong sign.

    dPI_dU[0] = Chi_i - 0.5 * Kappa_i * sq_vel_i;

  • The derivative of the star pressure retains the previous approximation by default, for compatibility (reported in Bug in HLLC Jacobian calculation #2198). USE_EXACT_HLLC_JACOBIANS= YES selects the exact derivative with wave-speed estimates held fixed, for ideal and general gases. The AUSM/SLAU numerical-Jacobian option remains separate. With the shared moving-wall correction from Fix the incompressible moving-grid Jacobian and several smaller mean-flow bugs #2945, both modes converge in the translating-airfoil checks; these results do not establish a general robustness ordering.

    dpStar_dU[iVar] = Density_i * (sR - ProjVelocity_j) * dSm_dU[iVar];

  • The ideal-gas supersonic Jacobians now use the reconstructed face enthalpy, matching the flux. The previous energy-based helper could give a different enthalpy after MUSCL reconstruction; both Jacobian modes shared this mismatch. General-gas HLLC already uses the enthalpy-based helper.

  • Verification access is private to the unit tests; no public wave-speed override is exposed. Configuration tests check that the HLLC and AUSM/SLAU Jacobian options are independent.

  • Unit tests compare the Jacobians with finite differences in 2D and 3D, on fixed and moving grids, and check independently reconstructed ideal-gas face pressure and enthalpy in both Jacobian modes.

The HLLC regression cases all use fixed grids. The test case is in the comments below.

This PR builds on #2942. It is part of the split of #2941.

Related Work

PR Checklist

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
  • My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
  • I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
  • I have added a test case that demonstrates my contribution, if necessary.
  • I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

🤖 Generated with Claude Code

https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6

- Moving grids: the wave speeds, the pressure work of the moving face in the supersonic branches and the term
  V_g (p* - p)/(s - S_M) of the star-state energy were missing or inconsistent, so HLLC was not Galilean
  invariant (NACA0012 M 0.8, fixed grid vs still air and a translating grid: CL -12 %). The flux is now
  computed in the frame of the face (u' = u - w n, H' = H - w u_n + w^2/2) and transformed back
  (F = B F', J = B J' T).
- Jacobians: sign of the velocity term of the general-gas dp/drho. The star-state Jacobians keep by default the
  density of the other state in dp*/dU (inexact, #2198), which was more robust with inexact linear solves in our
  tests; USE_ACCURATE_FLUX_JACOBIANS= YES now gives the exact derivative for fixed wave speeds.
- Unit test hllc_jacobian_tests.cpp: finite differences in 2D and 3D, fixed and moving faces, both sides of
  the contact, all branches, default and exact.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6
@rois1995

Copy link
Copy Markdown
Contributor Author

Test case

NACA0012, Mach 0.8, angle 1.25 deg, inviscid, HLLC with second-order reconstruction, 2 MPI ranks. Files and scripts: pr2.

  • Fixed grid (reference): air at Mach 0.8.
  • Moving grid: still air, and the grid moves at the same speed in the opposite direction.

The two must give the same flow. Before this PR (develop with PR 1), the moving case stops converging at a density residual of -4.2, with 12 % less lift (CL 0.2954 instead of 0.3351) and 2.4 % less drag. With this PR it reaches -11 in 5294 iterations, with CL 0.3359 and CD 0.02315, on top of the fixed-grid curves (CL 0.3351, CD 0.02318, -11 in 5148 iterations).

HLLC translating grid

The Jacobians are checked by the new unit test hllc_jacobian_tests.cpp.

@joshkellyjak joshkellyjak self-assigned this Oct 12, 2026

@joshkellyjak joshkellyjak left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/*--- Jacobians for the wave speeds held fixed (in the frame of the face). p* = rho_j (u_j - sR)(u_j - sM) + p_j = rho_i (u_i - sL)(u_i - sM) + p_i, so its derivative contains rho_j (sR - u_j) dsM/dU (or rho_i (sL - u_i) dsM/dU). By default the density of the other state is kept, as in the previous implementation (an inexact transcription of the reference, see #2198): in our tests it was more robust with inexact linear solves than the exact derivative, which USE_ACCURATE_FLUX_JACOBIANS= YES gives. ---*/

Can you provide the results of these tests to justify this?

Comment thread Common/src/CConfig.cpp Outdated
/*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the Lax-Friedrichs method. \ingroup Config*/
addDoubleOption("LAX_SENSOR_COEFF", Kappa_1st_Flow, 0.15);
/*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) \ingroup Config*/
/*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2), exact Jacobians (for fixed wave speeds) for HLLC \ingroup Config*/

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given the two behaviours here for AUSM/SLAU and HLLC, I think it would be better to split these into two separate options.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed. HLLC now has USE_EXACT_HLLC_JACOBIANS, defaulting to NO. USE_ACCURATE_FLUX_JACOBIANS retains its AUSM+up(2)/SLAU(2) behavior. I updated the constructors, documentation and tests, including all four combinations of the two options. “Exact” still means the wave-speed estimates are held fixed. The configuration and finite-difference checks pass in serial and on two MPI ranks.

* \param[in] active - Use the fixed values.
* \param[in] val_sL, val_sR - Left and right wave speeds.
*/
inline void SetFixedWaveSpeeds(bool active, su2double val_sL, su2double val_sR) {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If these are only used for the unit test they should not be in the public API

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed SetFixedWaveSpeeds from both public HLLC APIs. A friend test accessor sets the private verification state, so the finite-difference tests can still hold the wave speeds fixed.

const su2double rhoPStar_L = useAccurateJacobian ? Density_i : Density_j;



Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

trailing lines

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed the extra blank lines.

const su2double rhoPStar_L = useAccurateJacobian ? Density_i : Density_j;



Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed here too.

@rois1995

Copy link
Copy Markdown
Contributor Author

Merged the updated #2942 in d2701a8, including the latest develop merge and the under-relaxation review changes. I am repeating the HLLC Jacobian comparison on this updated code.

@rois1995 rois1995 changed the title Fix the HLLC scheme on moving grids and its Jacobians [WIP] Fix the HLLC scheme on moving grids and its Jacobians Oct 12, 2026
Reuse the enthalpy-based flux Jacobian helper in both ideal-gas supersonic branches. Check independently reconstructed face pressure and enthalpy in both Jacobian modes, 2D/3D, and on fixed/moving faces. Remove the general robustness claim from the code comments.
@rois1995

Copy link
Copy Markdown
Contributor Author

@joshkellyjak, I repeated the comparison on the current PR after the develop/parent updates and isolated a separate Jacobian bug before drawing conclusions.

Pushed the fix and regression test. In the ideal-gas HLLC supersonic branches, the flux uses reconstructed face enthalpy, but the old Jacobian helper recomputed enthalpy from energy and velocity. MUSCL reconstruction can make these values differ. Both exact and approximate modes share these branches. The correction reuses the existing enthalpy-based helper; general-gas HLLC already uses it.

The regression test covers independently reconstructed face pressure and enthalpy in both modes, 2D/3D, and on fixed/moving faces. Both HLLC test cases pass, with 784 assertions.

After the correction, the translating ideal-gas NACA0012 case (relative Mach 0.8, angle of attack 1.25 degrees, four MPI ranks) gives:

Jacobian CFL FGMRES settings Result
Approximate 5 LU-SGS, maximum 5 iterations, tolerance 1e-3 Converges at iteration 5296
Exact 5 LU-SGS, maximum 5 iterations, tolerance 1e-3 Produces NaNs after the last recorded iteration 1424
Exact 2 LU-SGS, maximum 5 iterations, tolerance 1e-3 Does not converge within 20,000 iterations; density residual remains near -2.28
Exact 5 ILU(0), maximum 50 iterations, tolerance 1e-6 Converges at iteration 5210

The convergence target is a log10 density residual of -11. All runs start from scratch with the same mesh and case configuration, apart from the Jacobian and linear-solver/CFL settings shown. "Exact" here means differentiation with the wave-speed estimates held fixed.

For comparison, before the enthalpy correction the approximate loose run converged at 5293, whereas the exact tight run produced NaNs after the last recorded iteration 1421. The correction therefore restores the tight-solve result while the difficulty with loose solves remains.

This supports the robustness observation for this particular case, rather than a general claim. I removed the broad claim from the code comments and qualified the PR description accordingly. These convergence results cover ideal gas only. The general-gas unit checks use ideal-gas thermodynamics through that implementation; an actual real-gas EOS/convergence comparison has not been run.

@rois1995

rois1995 commented Oct 12, 2026 •

Copy link
Copy Markdown
Contributor Author

I ran additional local checks against develop to assess the real-gas behavior and the moving-grid correction.

Scope clarification: the bump case is a Galilean-invariance check. The entire channel, including the bump and all boundaries, is translated, and the prescribed fluid velocities receive the same boost. The bump remains attached to its wall; this is the same flow viewed from two inertial frames. It is not intended as a simulation of a bump sliding along a stationary wall, or as the primary validation of moving-body startup.

The convergence case uses the existing nicf/edge/edge_PPR.cfg Peng–Robinson MDM configuration: Mach 1.7, approximately 15,000 mesh points, four MPI ranks, no multigrid, CFL 5, and a density residual target of 10^-11. I compared develop (18cc7a1e) with this PR (4960dd03), using both approximate and exact HLLC Jacobians. "Exact" here means differentiation with the wave-speed estimates held fixed.

On the fixed grid, develop converged in 1,794 iterations and both PR modes in 1,796. The converged fields agree closely, so this case shows no material fixed-grid improvement or regression.

Fixed-grid convergence

For the moving-grid comparison, I added the same 20 m/s translation to the grid, freestream, and inlet velocities, preserving the thermodynamic states and relative flow. The force normalization was also kept consistent. All three moving runs started from the same boosted, converged fixed-grid solution:

Version Iterations Lift change from fixed grid Drag change from fixed grid
develop 835 +1.747% −5.187%
PR, approximate Jacobian 553 −0.00123% −0.00098%
PR, exact Jacobian 553 −0.00123% −0.00098%

The PR substantially reduces the dependence on the reference frame. The remaining discrepancy is small but nonzero. Exact and approximate Jacobians give essentially identical results in this comparison.

Moving-grid convergence from identical boosted restart

In both plots, the wide translucent stroke highlights the approximate curve for visibility; its width has no uncertainty meaning.

There are important startup limitations:

  • From freestream with MUSCL, develop and both PR modes fail. Tightening the linear solve or reducing CFL to 1 did not resolve this.
  • At first order and CFL 5, develop converges in 751 iterations, while both PR modes fail.
  • The PR's approximate first-order run at CFL 1 remains finite, but does not converge within 20,000 iterations.

These results support improved frame consistency, but do not support a general robustness claim. Subsequent tracing identified the immediate startup failure in the shared characteristic outlet boundary: during the transient, a locally incoming acoustic characteristic activates the prescribed 1,000 Pa pressure and the boundary construction produces negative density. A standalone replay reproduces the invalid boundary state in both fixed and translating frames. Changing only the outlet pressure to 1.4 MPa permits both PR modes to converge from freestream, and the converged MUSCL solution matches the earlier moving-grid restart result.

The original startup failures should therefore be read as diagnostics of this extreme outlet-pressure setup, rather than evidence by themselves of an HLLC flux defect or a representative moving-body robustness limitation. The frame-invariance comparison remains useful. I am checking startup separately with a moving-body case that has no prescribed-pressure outlet.

Additional local unit checks exercise the production Peng–Robinson and van der Waals models in 2D/3D, on fixed and moving faces. The combined HLLC checks passed 916 assertions, including finite-difference checks of the exact Jacobians with wave speeds held fixed. These checks validate the tested flux derivatives, rather than full-solver convergence.

Experimental shared-wall corrections did not pass convergence checks and were reverted. No additional solver changes or regression testcase have been pushed for this investigation.

The case configurations, plots, and findings are published separately in my fork. Metadata records the tested revisions, mesh identity, and instructions for generating the common boosted restart. The raw restart state is not included.

@rois1995

rois1995 commented Oct 12, 2026 •

Copy link
Copy Markdown
Contributor Author

The boundary follow-ups are now separated from this HLLC change:

The results below describe the initial outlet checkpoint at 9307f8d; the current implementation and further validation are linked above.

For the characteristic outlet, the linear finite-amplitude state projection is replaced by reconstruction along the EOS isentrope, integrating density and normal velocity with the local sound speed. Expansion stops at the sonic trace when a prescribed pressure would require incoming information from a supersonic outlet. Already-supersonic outflow copies the interior state. Invalid thermodynamic states and backflow give explicit errors. The implicit convective boundary Jacobian differentiates this reconstruction.

Validation used two MPI ranks, first-order HLLC, implicit FGMRES/ILU and CFL 20. The straight duct starts at Mach 0.3 with outlet pressure 5% below its initial pressure. The smooth converging-diverging nozzle has throat/inlet area ratio 0.8, starts at Mach 0.9 and prescribes inlet density/velocity; its thermal inlet state comes from outgoing characteristics. These are fixed-wall physical flow checks, separate from the translating-bump Galilean-invariance diagnostic.

Case develop repaired outlet
Air, subsonic duct 1571 1403
MDM/Peng–Robinson, subsonic duct 1595 1467
Air, nozzle 1062 1062
MDM/Peng–Robinson, nozzle NaN at startup 3616

Numbers are iterations to log10 density residual −10. The repaired duct outlets remain subsonic and recover the requested pressure (relative errors about 1e−11 for air and 1e−13 for MDM). In the nozzles, the configured 1000 Pa drives startup; final outlet flow is supersonic (Mach about 1.59 for air and 1.42 for MDM), so that pressure is no longer imposed.

Configurations, nozzle mesh, portable runner and summary are available. The duct convergence plot compares the baseline and repaired outlet. The baseline was develop 056199c; the initial repaired binary was rerun on all four selected cases.

The focused boundary tests pass 1364 assertions, covering ideal gas, Peng–Robinson and van der Waals, 2D/3D, moving-grid boosts, sonic transition, supersonic identity, frozen turbulent kinetic energy, and convective-flux Jacobian finite differences. All three Meson test suites pass. Test source.

Limitations: this is tested single-phase outflow, not a general backflow boundary. A total-pressure/temperature-inlet MDM nozzle variant stays finite but stalls near density residual −8.8; a straight duct approaching sonic flow also converges slowly. Those variants remain in the evidence summary, so I am not claiming universal startup robustness or faster convergence. Further mesh, AD and OpenMP checks are now documented in #2975. Full adjoint flow runs remain unchecked, and the existing viscous Jacobian approximation is unchanged. The reservoir-inlet residual floor was subsequently isolated to EOS inversion precision, as described in the separate follow-up linked above.

@rois1995

Copy link
Copy Markdown
Contributor Author

The shared compressible wall correction is now pushed to #2945 in e0d3ee8d72. Its dedicated boundary/multigrid frame-covariance regression and the #2945 unit suites pass. The earlier local combined runs converge with the original real-gas bump outlet setting in both Jacobian modes; the ideal-gas bump and translating NACA exact/loose-linear-solve run also converge. This is a separate wall-boundary correction, not another HLLC flux change.

@rois1995

Copy link
Copy Markdown
Contributor Author

The explicit supersonic-outlet setup fix is now in #2974. Its boundary regression reproduces the original general-gas HLLC segfault and passes 8,580 assertions after supplying secondary EOS data in both convective and viscous calls. The checks cover ideal/Peng–Robinson, 2D/3D, fixed/moving grids, and Euler/Navier–Stokes. This is separate from the characteristic static-pressure outlet used by the bump; its nonlinear repair remains on a separate branch.

@rois1995

Copy link
Copy Markdown
Contributor Author

@joshkellyjak, I removed the broad robustness statement. Subsequent tracing found a shared moving-wall projection bug: with the correction from #2945, the exact-Jacobian translating NACA run also converges with the loose linear solve, reaching the target at iteration 5155.

The real-gas fixed-grid comparison gives 1794 iterations on develop and 1796 in both PR modes. The boosted-restart comparison also gives matching results between the two PR modes. The separate wall and outlet findings, tests and limitations are documented in the follow-up comments.

I am therefore retaining the approximate default for compatibility, without claiming that it is generally more robust than the exact fixed-wave-speed derivative.

@rois1995

Copy link
Copy Markdown
Contributor Author

Propagated #2942’s AD density lifetime fix and its solver regression in 8b28ba8eab. Both previously failing file comparisons now pass without reference changes; the forward-AD suite and focused reverse-AD/ASAN regression pass. Details.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants