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 crashed. It is now limited to 10, with a warning.

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); with it, it warns and uses 10 levels.

[joshkellyjak]

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

[joshkellyjak]

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

[joshkellyjak]

const auto

[rois1995]

Merged the updated #2944 in 5ca9969 to propagate the latest #2942/develop changes through the stack.

[rois1995]

Propagated the ideal-gas HLLC supersonic face-enthalpy Jacobian fix and the regression test covering both Jacobian modes from #2943 (4960dd03be). The HLLC unit tests pass with 784 assertions. The convergence comparison and its scope are documented in #2943.