Fix the internal energy used by the under-relaxation factor

[rois1995]

Proposed Changes

The implicit compressible solver limits the update of each point with an under-relaxation factor. To decide how much to limit, it compares the internal energy before and after the update. This internal energy was computed wrongly: the kinetic energy was not divided by the density.

SU2/SU2_CFD/src/solvers/CEulerSolver.cpp

Lines 2555 to 2556 in a6b7496

	e_old -= 0.5 * pow(nodes->GetSolution(iPoint, jVar), 2);
	e_new -= 0.5 * pow(nodes->GetSolution(iPoint, jVar) + LinSysSol(iPoint, jVar), 2);

In cases run with dimensional values the wrong internal energy can become negative. When it does, the point gets no update at all and stays frozen. In an air nozzle about a quarter of the points were frozen, so the case could not converge.

The fix is 5 lines: the kinetic energy is divided by the density, using the old and the new density.

Almost all compressible implicit regression cases converge a little differently, so most regression values change. They are updated from the CI.

This is the first of four PRs that replace #2941, to make the review easier. The other three build on this one:

• Fix the HLLC scheme on moving grids and its Jacobians #2943: HLLC on moving grids;
• Fix L2Roe, LMRoe and the general-gas Roe scheme in 3D #2944: Roe-type schemes in 3D;
• Fix the incompressible moving-grid Jacobian and several smaller mean-flow bugs #2945: smaller fixes.

Related Work

• Detect non-physical solutions in a statistic way instead of just P or T <= 0 #2830 added the under-relaxation factor.
• Fix mean-flow solver bugs: under-relaxation, HLLC on moving grids, moving-grid Jacobians and others #2941 (replaced by this PR and the three that follow).

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

[rois1995]

Test cases

Both cases: axisymmetric air nozzle, dimensional, 2 MPI ranks, develop against this PR. Files and scripts: pr1 (see the README).

Restart from a converged solution (SST). The start solution was converged with this PR. Develop does not stay converged: the density residual rises from -12 to -10.2 in 300 iterations, because the points with a "negative internal energy" are frozen. With this PR the residual stays at -12 and the run stops after 10 iterations.

Start from rest (SA, CFL adaptation). Develop stays at a density residual of -2 for 20 000 iterations, with a wrong outlet pressure (5.1 bar). With this PR it reaches -11 in 11 647 iterations (outlet pressure 1.64 bar).

[bigfooted]

@rois1995 if this PR is finished can you update the regression tests? looks like it's a small PR that can be merged soon.

[rois1995]

Pushed 3062a04822 to materialize rho_old and rho_new as const su2double. The previous const auto declarations retained CoDi expressions referencing expired temporary values from GetSolution(), causing the reported stack-use-after-scope.

Added an actual-solver AD regression checking relaxation factors and density derivatives for density- and energy-limited updates. The forward-AD unit suite passes, and the new under-relaxation regression passes in reverse AD with AddressSanitizer. Both previously failing direct-differentiation file comparisons now match the existing references with zero tolerance; no reference files were changed.

The fix and test are also propagated through #2943, #2944 and #2945. Full CI will confirm the remaining configurations.

[bigfooted]

There are some concerns regarding the EPS floor. I am not too worried about it and will approve the PR, just something to test later maybe and see if this is better or not:

Why the floor blocks the update. Say the full update would push rho_new to or below 0.

The density ratio |Δrho|/rho_old is then at least 1. With MAX_UPDATE_FLOW below 1, the density limit alone gives a factor of allowable/ratio, which is small but nonzero.
The energy term is evaluated at rho_new = EPS. Its kinetic part, |rho·u|²/(2·EPS), is huge, so the energy ratio is astronomically large.
The combined factor drops below 1e-10 and is zeroed. The point gets no update, which is the same frozen-point symptom the PR is fixing.
So a legitimate "take a smaller step" case becomes "don't move at all". The limiter should shorten the step, not cancel it.

Options, best first

Linearize the energy change. This is what I'd do. With e = E − m²/(2ρ), the first-order change is:

Δe ≈ ΔE − u·Δm + ½|u|²·Δρ

It needs no division by rho_new and no floor, and it stays bounded.
The limiter only has to size the step, so first-order accuracy is enough. A large update, where the second-order terms matter, is exactly the case being limited.
It is cheaper and has cleaner AD behaviour.
Caveat: it can underestimate the change when the density drops sharply. The density ratio already covers that case.

It's a small code change but a larger validation job.

Code. It is a rewrite of the loop body in CEulerSolver::ComputeUnderRelaxationFactor, about 10 lines.

const su2double rho = max(nodes->GetSolution(iPoint, 0), EPS);
su2double e = nodes->GetSolution(iPoint, nVar - 1);   // rho*e
su2double de = LinSysSol(iPoint, nVar - 1);
for (auto d = 1u; d <= nDim; ++d) {
  const su2double u = nodes->GetSolution(iPoint, d) / rho;
  e -= 0.5 * rho * u * u;
  de -= u * LinSysSol(iPoint, d) - 0.5 * u * u * LinSysSol(iPoint, 0);
}
ratio = fmax(fabs(LinSysSol(iPoint, 0)) / rho, fabs(de) / max(e, EPS));

Energy change. de is the Δe formula from before. Only rho is floored, and it is a divisor of the old state, so it can't be driven to EPS by the update.
Old energy. e is the same rho*e the PR computes.
Unchanged. The rest of the function stays as it is: the fmin(allowable/ratio, 1) line, the 1e-10 cancel and SetUnderRelaxation.
What changes numerically.

For small updates it matches the PR's formula to second order, so well-behaved cases barely move.
For updates where rho_new is far from rho_old, it differs from the exact e_new − e_old. It no longer produces the huge ratio that freezes the point.
Cases that currently freeze on a near-vacuum rho_new start moving again, scaled by the density limit.
Work around the code.

Unit test. Extend under_relaxation_ad_tests.cpp. The expected factor and derivative change, and a new case should drive rho_new to 0 and assert that the factor is positive and equals the density-limited value. The AD/DD derivative checks need new analytic values, so recompute them.
AD. The expression uses only active values and has no data-dependent branch beyond the max calls, so it should record fine. The test is what confirms it.
Regression references. Every implicit compressible case shifts a little, as it does in this PR. Regenerate on top of this PR's references so the diff is attributable.
Convergence check. Run a few known-sensitive cases, such as the dimensional nozzle from this PR, a strong-shock case and a startup from poor initial conditions. Compare iteration counts to the exact-formula version, to confirm that the underestimate in the density-drop case doesn't allow overshoots.
Docs. MAX_UPDATE_FLOW should state that the energy limit uses a first-order estimate.
Risk. The one real risk is the underestimate. If the linearized de is too small when rho drops and |u| is large, the limiter allows a step that the exact formula would have shortened, which could cause overshoot in strong rarefactions. The density ratio covers most of this, so I expect it to be fine, but step 4 is what settles it.

[bigfooted]

Looks good, see my comment. Maybe test later and see if it helps any cases to do a linearization.

[rois1995]

@bigfooted, thanks for the comment! Should we test that now before merging?

[bigfooted]

@rois1995 sure you can test now and see if it helps. if it is a very clear win you can add it now, if it is inconclusive we can merge this now.

[rois1995]

@bigfooted I tested the linearized version against this PR.

Smooth cases: same result.

Case	This PR	Linearized
Air nozzle (axisymmetric, dimensional), restart, SST, CFL 1000	392 it	392 it
Air nozzle from rest, SA, adaptive CFL	8280 it	10957 it, same outlet pressure
RAE2822 SA, JST, multigrid, CFL 1000	2200 it	2241 it, same CL/CD
Blunt body M6, 1st order, CFL 6, freestream start	1742 it	1767 it, same CD
Ramp M2, freestream start, adaptive CFL and fixed CFL 50/200	149/169/123 it	151/170/124 it

In these cases no point is frozen with this PR. On the nozzle from rest, the two histories are the same up to about iteration 5000. After that, only the slow end of the convergence differs.

Hard start-ups: the linearized version is worse. Blunt body M6 from freestream, 6000 iterations (final log rms[Rho], iterations with non-physical points in brackets):

Run	This PR	Linearized
1st order CFL 10	converged (1297 it)	-0.53 (3)
1st order CFL 20	0.69	diverges (4040)
1st order CFL 50	0.80	2.29 (501)
1st order CFL 100	diverges (67)	2.48 (104)
1st order CFL 200	2.63	2.94 (1110)
MUSCL CFL 5	converged (4627 it)	converged (4430 it)
MUSCL CFL 10	-0.99	-1.00 (7)
MUSCL CFL 20	0.21	diverges (1800)
MUSCL CFL 50	0.48	diverges (1980)

This PR: 2 converged, 6 stalled, 1 diverged. Linearized: 1 converged, 5 stalled, 3 diverged, and non-physical points in 8 of 9 runs.

You were right about the freeze. With this PR, a point freezes only when the full update makes the density zero or negative. This happened only in these hard runs. At 1st order CFL 10 the run still converges, and the linearized version stalls there.

The reason: the exact change of rho*e is the linear one minus |d(rho u) − u d(rho)|² / (2 rho_new). This term always lowers the internal energy. The linear estimate misses it, so large velocity updates can make the energy negative. The density limit does not catch this.

The linearized version would also change 7 of the 13 serial regression cases that I ran by hand.

So I would merge this PR as it is.