add support for axisymmetric jets with SA turbulence model

[bigfooted]

Add support for axisymmetric jets with SA turbulence model. SA needs an additional submodel:
https://www.researchgate.net/publication/349801309_A_New_l2_Term_for_the_Spalart-Allmaras_Turbulence_Model_Active_in_Axisymmetric_Flows

[Syphonicc]

I had a look at this. Reporting what I found and what I have so far, rather than opening a PR blind the issue has no description, so please correct me if this is not what you meant.

What is already there

SA with AXISYMMETRIC= YES is not blocked: the guard "Axisymmetry is currently only supported for KIND_TURB_MODEL chosen as SST" was removed in d90e8e1 (Jan 2024), which added ResidualAxisymmetricDiffusion to CSourceBase_TurbSA. The axisymmetric nozzle case runs with KIND_TURB_MODEL= SA without error.

What looks incomplete

SA adds only the diffusion part of the axisymmetric contribution. The two comparable transported scalars add convection as well:

scalar	axisym. convection	axisym. diffusion	Jacobian
SST k, omega (turb_sources.hpp)	-yinv*V*rho*v*phi	yes	yes
species Y (species_sources.cpp)	-yinv*V*rho*v*phi	yes	yes
SA nu_tilde	missing	yes	none

SU2 discretizes the SA convective term in divergence form (Flux = a0*nue_i + a1*nue_j, no density). In cylindrical coordinates that divergence is d(u*nue)/dx + (1/r) d(r*v*nue)/dr, so the 2D planar fluxes miss (1/r)*v*nue. Equivalently: for a uniform nue in an incompressible axisymmetric flow the convective residual should vanish, because du/dx + dv/dr + v/r = 0; without the term it evaluates to -nue*v/r per unit volume. That error is proportional to the radial velocity, so it is small in a boundary layer and not small in a jet or a diverging nozzle.

Change I have on a branch

In CSourceBase_TurbSA, mirroring the SST naming:

• ResidualAxisymmetricDiffusion -> ResidualAxisymmetricConvectionDiffusion, adding Residual[0] -= yinv*Volume*v*nue and Jacobian_i[0][0] -= yinv*Volume*v (velocity only, no density, to match SA's density-free convective flux).
• The call is moved to after Jacobian_i[0][0] *= Volume;, otherwise the new Jacobian term would be scaled by the volume twice. Residual behaviour of the existing diffusion term is unchanged by the move.
• Coord_i[1] and V_i[idx.Velocity()+1] are registered as preaccumulation inputs when axisymmetric. Coord_i[1] was already used by the existing diffusion term without being registered, so the discrete adjoint of axisymmetric SA was missing that dependency.

Verification so far

• Non-axisymmetric is untouched: everything new sits behind if (axisymmetric), and turb_naca0012_sa at iteration 5 still gives -12.037372, -16.384159, 1.080346, 0.018385, 20, -3.456391, 20, -4.641246, 0 against the committed -12.037323, -16.384159, 1.080346, 0.018385, 20, -3.457510, 20, -4.641262, 0 (serial GCC 15.2 build, differences only in the linear-solver columns).
• The term is active in the axisymmetric nozzle case (axisymmetric_rans/air_nozzle switched to SA, cold start, 1000 iterations): rms[nu_tilde] -7.327 before, -7.155 after.

That second point only shows the term does something, not that it is right. Before I open a PR I would like to validate it against a 3D wedge of the same geometry, which is the only check that really settles the sign and the factor, and add an axisymmetric SA regression case so this path is covered in CI. Two questions:

1. Is this the gap you had in mind, or were you after a validation/tutorial jet case?
2. For the regression entry, would you prefer an SA variant of the existing air_nozzle case (mesh already in the TestCases repo), or a dedicated jet case (new mesh)?

[bigfooted]

SA cannot simulate axisymmetric jets, it needs a separate submodel:
https://www.researchgate.net/publication/349801309_A_New_l2_Term_for_the_Spalart-Allmaras_Turbulence_Model_Active_in_Axisymmetric_Flows

[Syphonicc]

Thanks, that is the piece I was missing I had been looking at the discretization rather than the model. For anyone reading later, that paper is documented on the NASA TMR as SA-TC (transverse curvature free-shear correction), Spalart & Garbaruk, Flow, Turbulence and Combustion 107:245-256 (2021).

The term added to the RHS of the SA equation is

max(1 - r, 0)^3 * c_b3 * lambda_2 * nu_hat / sigma,    c_b3 = 6, sigma = 2/3

with r the same as in fw, so it switches off near walls and is fully active in free shear, and lambda_2 the middle eigenvalue of the Hessian of nu_hat. The mechanism is that for a plane jet the out-of-plane Hessian component is zero and the term vanishes, while for a round jet the azimuthal component of the Hessian of an axisymmetric field is nu_hat_r / r, which is not - that is what fixes the round-jet/plane-jet anomaly.

I would like to implement it. From reading the code, the pieces are:

1. A new SA_OPTIONS= TRANSVERSE_CURVATURE entry and the corresponding flag in SAParsedOptions.
2. The Hessian of nu_tilde, which does not exist anywhere in SU2 today. The route that does not require new infrastructure is the existing auxiliary-variable gradient: set nAuxVar = nDim for the SA variable, store the components of grad(nu_tilde) in AuxVar, and call SetAuxVar_Gradient_GG/LS, which already has an MPI communication type. The result is the Hessian at each point.
3. lambda_2 itself: in 2D planar, the eigenvalues of the 2x2 in-plane Hessian plus a zero third one; in 2D axisymmetric, the 2x2 block plus nu_tilde_r / r as the out-of-plane eigenvalue; in 3D the analytic eigenvalues of a symmetric 3x3. Middle = sum - min - max in all cases.
4. The source term and its Jacobian, which is straightforward since the term is linear in nu_tilde.

One thing I would like your opinion on before I get too far. Neither the paper nor the TMR page prescribes how to discretize the second derivatives, and a Green-Gauss gradient of a Green-Gauss gradient is noisy on unstructured meshes, which is exactly where lambda_2 would be evaluated in a shear layer. If you know how this was done in the reference implementation, or have a preference (Green-Gauss twice, least-squares on the gradients, something else), that would save me guessing.

For validation I would target the round jet spreading rate against the plane jet, which is what the paper is about. That needs a jet mesh, which means something in the TestCases repo; tell me if you would rather I propose a case there first.

---

Separately, and not part of this: while reading CSourceBase_TurbSA I noticed the axisymmetric contribution only includes the diffusion part. SST and the species solver both also add the convection term -(1/y) * rho * v * phi, and SU2 writes the SA convective flux in divergence form, so the cylindrical divergence appears to be missing (1/y) * v * nu_tilde. Equivalently, for a uniform nu_tilde in an incompressible axisymmetric flow the convective residual should vanish because du/dx + dv/dr + v/r = 0, and without that term it does not. I have a patch for it but have not validated it against a 3D wedge yet, so I am not claiming it is right. Should I open a separate issue for that, or is it deliberate?