Fix the incompressible moving-grid Jacobian and several smaller mean-flow bugs

[rois1995]

Proposed Changes

Each fix is in its own commit.

• Incompressible solver on moving grids. The Jacobian term of the grid motion was added once per variable instead of once. The results of each step do not change, but the inner iterations of unsteady cases converge much worse. In a pitching airfoil, the lift reached values of several hundred because the time steps did not converge. The term is now added once, including the effect of a variable density.

  SU2/SU2_CFD/src/numerics/flow/convection/fds.cpp

  Lines 252 to 265 in a6b7496

  	for (iVar = 0; iVar < nVar; iVar++) {
  	Flux[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]);
  	/*--- Jacobian contributions ---*/
  	/*--- Implicit terms ---*/
  	if (implicit) {
  	for (iDim = 0; iDim < nDim; iDim++){
  	Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i;
  	Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j;
  	}
  	Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i;
  	Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j;
  	}
  	}

• Source terms that cannot be combined. Some source terms share one internal object, so when two are active, one silently replaces the other. For example, body force + gravity gave twice the body force and no gravity. The first source slot is one object for all of them:

  SU2/SU2_CFD/src/drivers/CDriver.cpp

  Lines 1843 to 1857 in a6b7496

  	numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncBodyForce(nDim, nVar_Flow, config);
  	else
  	numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceBodyForce(nDim, nVar_Flow, config);
  	}
  	else if (incompressible && (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE)) {
  	numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncStreamwise_Periodic(nDim, nVar_Flow, config);
  	}
  	else if (incompressible && (config->GetKind_DensityModel() == INC_DENSITYMODEL::BOUSSINESQ)) {
  	numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceBoussinesq(nDim, nVar_Flow, config);
  	}
  	else if (config->GetRotating_Frame() == YES) {
  	if (incompressible)
  	numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncRotatingFrame_Flow(nDim, nVar_Flow, config);
  	else
  	numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config);

  These combinations now stop with a clear error.

• Turbomachinery Riemann boundary. It read the settings of the Giles boundary instead of its own, so it crashed when no Giles marker was present.

  SU2/SU2_CFD/src/solvers/CEulerSolver.cpp

  Line 5800 in a6b7496

  switch (config->GetKind_Data_Giles(Marker_Tag)){

  It now reads its own settings.

• Reconstructed face states (MUSCL). Pressure and density are now checked to be positive before the fluid model uses them. Before, the fluid model was called first:

  SU2/SU2_CFD/src/solvers/CEulerSolver.cpp

  Lines 1952 to 1968 in a6b7496

  	if (!ideal_gas || low_mach_corr) {
  	ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i);
  	ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j);
  	}
  	/*--- Low-Mach number correction. ---*/
  	if (low_mach_corr) {
  	LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j);
  	}
  	/*--- Check for non-physical solutions after reconstruction. If found, use the
  	cell-average value of the solution. This is a locally 1st order approximation,
  	which is typically only active during the start-up of a calculation. ---*/
  	bool neg_pres_or_rho_i = (Primitive_i[prim_idx.Pressure()] < 0.0) || (Primitive_i[prim_idx.Density()] < 0.0);
  	bool neg_pres_or_rho_j = (Primitive_j[prim_idx.Pressure()] < 0.0) || (Primitive_j[prim_idx.Density()] < 0.0);

• Incompressible minimum pressure. With MPI it took the largest of the local minima.

  SU2/SU2_CFD/src/solvers/CIncEulerSolver.cpp

  Line 2238 in a6b7496

  SU2_MPI::Allreduce(&minP, &MinP, 1, MPI_DOUBLE, MPI_MAX, SU2_MPI::GetComm());

  It now takes the smallest.

• Axisymmetric source. One Jacobian term was always zero, because 1/2 is an integer division in C++.

  SU2/SU2_CFD/src/numerics/flow/flow_sources.cpp

  Line 141 in a6b7496

  jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 1/2*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) );

• Green-Gauss gradient. A fixed-size array had no size check (mixtures with many species could overflow it).

  SU2/SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp

  Line 72 in a6b7496

  static constexpr size_t MAXNVAR = 20;

  It is now as large as the one in the limiters, with an error above it.

• MGLEVEL above 10 could exceed the fixed multigrid arrays. Configuration preprocessing now rejects these settings with the supported range. Unit tests preserve levels 0 and 10 and verify the diagnostic for 11.

• Compressible moving slip walls. The density/momentum/energy projection is now applied consistently to the state, residual, Jacobian and multigrid paths. A boundary/multigrid frame-covariance regression covers the correction. The associated convergence checks and separate outlet work are linked from [WIP] Fix the HLLC scheme on moving grids and its Jacobians #2943.

The test cases are in the comments below.

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

Related Work

• Multizone adjoints for turbomachinery #2446 (Turbomachinery Riemann boundary).
• Fix mean-flow solver bugs: under-relaxation, HLLC on moving grids, moving-grid Jacobians and others #2941 (split).

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

All runs: the build before this PR (PR 3) against this PR, 2 MPI ranks. Files and scripts: pr4.

Incompressible pitching airfoil. NACA0012, pitching by 5 deg, 60 inner iterations per time step. Before this PR the inner iterations stop converging (worst end-of-step pressure residual -2.05) and the lift reaches -178 and +388. With this PR the worst residual is -3.85 and the lift follows the motion.

Other fixes:

• Source terms: inviscid NACA0012 with a body force of (0, -50, 0) and gravity. Before this PR it gives exactly the same history as a body force of (0, -100, 0) alone, so gravity is lost. With this PR the run stops with an error that explains which options cannot be combined.
• Turbomachinery Riemann boundary: the 2D transonic stator with Riemann boundaries instead of Giles. Before this PR it crashes (segmentation fault); with it, it runs.
• MGLEVEL: turbulent flat plate with MGLEVEL= 12. Before this PR it crashes (stack smashing). The current implementation rejects unsupported levels during configuration preprocessing with a clear error. Unit tests check that levels 0 and 10 are retained and that 11 is rejected.

[joshkellyjak]

If this is in preprocessing I think it is better to throw an error rather than clamp with a warning to the screen

[rois1995]

Agreed. MGLEVEL > 10 now fails during configuration preprocessing with a diagnostic giving the supported range, instead of warning and clamping. Tests verify that levels 0 and 10 are retained and that 11 terminates with the diagnostic. All three configured test suites pass.

[joshkellyjak]

can you update the variable definition for variable_density so they agree in this file

[rois1995]

Both centered-flow constructors now use GetVariable_Density_Model() for variable_density.

[joshkellyjak]

const auto

[rois1995]

Changed both dRhodh declarations to const auto.

[rois1995]

The additional fixes are now in place:

• Compressible moving slip walls: the wall projection is applied consistently to the solution, residual, Jacobian and multigrid operations. The tests check that changing the reference frame gives equivalent boundary behaviour. The moving-grid convergence findings and their limitations are discussed in #2943.
• HLLC supersonic Jacobian: this PR includes the correction from [WIP] Fix the HLLC scheme on moving grids and its Jacobians #2943 to the derivative of the face-enthalpy term. The regression covers both the exact and approximate Jacobian modes. Explanation and tests.
• AD density lifetime: this PR also includes the fix from Fix the internal energy used by the under-relaxation factor #2942 that keeps the density intermediate alive during differentiation. The two previously failing file-comparison tests now pass with their original reference files. Forward-mode tests and the focused reverse-mode test with AddressSanitizer pass. Details.

The branch includes the updated parent PRs and develop changes. The wall and HLLC unit checks pass, and CI is currently passing on this PR's latest commit.

The pressure-outlet reconstruction is handled separately in #2975; it is not required by the wall correction here.