Fix the second-order convection of turbulence and species - #2956

Test cases

Files (configurations, histories, plot script) for each case: scalarMuscl. Meshes are from the TestCases repository. 2 MPI ranks, "develop" = a6b7496b74.

1. Values outside the allowed range: rans/naca0012/turb_NACA0012_sst.cfg, 10°, from the free stream, MUSCL_TURB= YES.

2. Bounded scheme with density ≠ 1: incomp_rans/rough_flatplate/rough_flatplate_incomp.cfg without roughness, BOUNDED_SCALAR, water density (998.2) with the viscosity scaled to keep the Reynolds number. It must match the run with density 1.

3. Limiter lost after LIMITER_ITER: rans/naca0012/turb_NACA0012_sa.cfg, restart, turbulence limiter VENKATAKRISHNAN_WANG, LIMITER_ITER= 50.

4. VAN_ALBADA_EDGE for turbulence: rans/flatplate/turb_SST_flatplate.cfg, MUSCL_TURB= YES. Develop stops with "Unknown limiter type"; this PR converges (CD 0.0028093).

5. Flow VAN_ALBADA_EDGE on the turbulence velocity: case 3 with SLOPE_LIMITER_FLOW= VAN_ALBADA_EDGE. CL 1.0876936 on develop, 1.0878103 with this PR.

6. Sliding interfaces: sliding_interface/bars_SST_2D with BOUNDED_SCALAR and KIND_INTERPOLATION= WEIGHTED_AVERAGE. With this PR the bounded scheme follows SCALAR_UPWIND: rms[Rho] in zone 1 after 500 iterations is -4.399, against -4.397 for upwind and -4.207 before.

7. Reconstruction gradient without MUSCL: case 4 with first-order turbulence and NUM_METHOD_GRAD_RECON= LEAST_SQUARES, 1 rank. The histories are identical; 0.238 s per iteration instead of 0.296 s.

New regression test turb_naca0012_sst_muscl_edge: develop stops with "Unknown limiter type". The test values change if fix 3 or fix 5 is undone.

can you check that the changed options are used in regression tests? I don't think we use van_albada_edge for instance for any species case or bounded scalar scheme with varying density?

@bigfooted added the missing species VAN_ALBADA_EDGE regression in 5aeb5f4, using variable density and BOUNDED_SCALAR in the venturi mixing case. I also added a stored turbulence limiter case that runs beyond LIMITER_ITER=3. The existing species2_primitiveVenturi_mixingmodel_boundedscalar already covers variable-density bounded transport. All four focused cases pass with two MPI ranks; the 44 unit tests pass too.

@rois1995 if you can fix the issues in the regression tests then we can finalize this.

@bigfooted There's still a problem with eigen that is outside the scope of this PR.

@bigfooted There's still a problem with eigen that is outside the scope of this PR.

that's probably the server being overloaded, it tries to get external packages but fails due to time-out:

fatal: remote error: GitLab is currently unable to handle this request due to load (ID a49f7cf55aa93bb9-LAX).

@rois1995 can you check this:

1. The scalar VAN_ALBADA_EDGE is not the flow's van Albada, and for turbulence it barely limits anything.
   The scalar kernel calls musclMusclLimiter::EDGE in SU2_CFD/include/numerics/util.hpp:

lim = (Δ² + 1e-6 + p·Δ) / (p² + Δ² + 1e-6)
The flow solvers use LimiterHelpers::vanAlbadaFunction(p, Δ, EPS), with EPS = 1e-16:

lim = Δ(p + Δ) / (p² + Δ² + 1e-16)
These differ in two ways:

The fixed 1e-6 also sits in the numerator, so as Δ→0 the limiter goes to 1 instead of 0. At an extremum the slope is left unlimited.
1e-6 is an absolute value. Non-dimensional ν̃ or k (e.g. with FREESTREAM_PRESS_EQ_ONE) is around 1e-8 to 1e-5, so Δ² ≪ 1e-6 and lim ≈ 1. The new turb_naca0012_sst_muscl_edge test therefore runs k with almost no limiting.
The PR's third fix applies this same kernel to the flow velocity that carries the turbulence. Δu is about 1e-3, so Δ² is about the size of the epsilon. The scalar face velocity is then still not the flow solver's own face velocity, which was the point of that fix.
Fix: call vanAlbadaFunction with EPS in the EDGE branch, or use an epsilon scaled to the variable. This changes the new test's reference values.

My Claude analysis found something:
The fix is to build the boundary mass flux with the upwind density, the same way the face itself is upwinded. Then the volume flux is u·n on both inflow and outflow, and the conservative flux takes the ghost density on inflow.

What's wrong now. BoundedScalarBCFlux always builds the mass flux from the interior density:


const su2double edgeMassFlux = density * DotProduct(nDim, velocity, normal); // ρ_i · (u_ghost·n)
This value is returned and passed to ComputeFlux, which turns it into volume fluxes by dividing by the density on each side:


a0 = fmax(0.0, massFlux) / rho.i; // outflow: ρ_i u·n / ρ_i = u·n ✓
a1 = fmin(0.0, massFlux) / rho.j; // inflow: ρ_i u·n / ρ_ghost = (ρ_i/ρ_ghost) u·n ✗
Inflow with a non-conservative scalar: the face carries (ρ_i/ρ_ghost)·u·n instead of u·n.
Inflow with a conservative scalar: the face carries ρ_i·u·n instead of ρ_ghost·u·n.
The PR's new BoundedScalarDivergenceFlux reproduces the same error, because it takes this mass flux as its input.

The fix, in CScalarSolver.hpp:


inline su2double BoundedScalarBCFlux(unsigned long iPoint, bool implicit, const su2double& density,
const su2double* velocity, const su2double* normal,
const su2double* densityGhost = nullptr, su2double weight = 1.0) {
const su2double rhoGhost = densityGhost ? densityGhost : density;
const su2double volumeFlux = GeometryToolbox::DotProduct(nDim, velocity, normal);
/--- Upwind density: the interior on outflow, the ghost on inflow. ---*/
const su2double edgeMassFlux = volumeFlux * (volumeFlux > 0.0 ? density : rhoGhost);
const su2double q = weight * BoundedScalarDivergenceFlux(edgeMassFlux, density, rhoGhost);
LinSysRes.AddBlock(iPoint, nodes->GetSolution(iPoint), -q);
if (implicit) Jacobian.AddVal2Diag(iPoint, -q / (Conservative ? density : 1.0));
return edgeMassFlux;
}

Gradient_Aux is the second gradient array a scalar solver keeps for MUSCL reconstruction. The PR stopped computing it when the scalar doesn't reconstruct, but each scalar variable still allocates it.

The fix. Make the variable use the same condition that CommonPreprocessing now uses. In the constructor, the generic config->GetMUSCL() is not set yet (it's a per-solver global parameter set during the iteration), so each derived class passes its own flag:

// CScalarVariable.hpp / .cpp
CScalarVariable(unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config, bool muscl);

CScalarVariable::CScalarVariable(..., const CConfig* config, bool muscl)
: CVariable(npoint, ndim, nvar, config),
Gradient_Reconstruction(muscl && config->GetReconstructionGradientRequired() ? Gradient_Aux : Gradient) {
Gradient.resize(nPoint, nVar, nDim, 0.0);
if (muscl && config->GetReconstructionGradientRequired()) {
Gradient_Aux.resize(nPoint, nVar, nDim, 0.0);
}

@bigfooted addressed the October 14 findings in eef97f4:

• Scalar/SIMD VAN_ALBADA_EDGE now calls the flow's vanAlbadaFunction with EPS. Added checks for extrema and small jumps; refreshed the SST and species edge-test references.
• Bounded boundary transport now upwinds density: interior on outflow, ghost on inflow. All callers pass the ghost density, and sliding donors retain their weights. Tests cover conservative/non-conservative corrections, explicit/implicit updates, zero flux and donor weights.
• Scalar variables allocate the auxiliary reconstruction gradient only when their own turbulence/species/heat MUSCL flag requires it. Poisson variables use the main gradient. Tests cover the constructor flags and shared/separate gradients.
• Clarified that only stored point limiters freeze after LIMITER_ITER; edge limiters are recomputed on every edge. The requested co-author attribution was already added in 5aeb5f4158.

Validation: 50 supported normal unit cases pass (75,383 assertions); the three new scalar tests pass on two MPI ranks and with two OpenMP threads (380 assertions per rank/run). All four focused MPI regressions pass. Formatting hooks pass. The local build disables MLPCpp, so its existing ANN fluid unit was excluded; CI will run the complete configured matrix.

The separate epsilon improvement is now open as #2978, with the priority label. It includes scaling, near-zero-reference, derivative and SA/SST regression tests, and targets this branch until #2956 merges.

@bigfooted the remaining eight x86_64 V&V/adjoint failures were reproduced with the binaries from CI run 38595536164. The shared edge helper also serves SIMD flow reconstruction, so its van Albada correction changes those cases. Updated their references in 48e871b; all eight pass with those same CI binaries, without changing tolerances. New unit files also carry author attribution (a3968c9).

The variable-relative scalar limiter epsilon follow-up is now #2978, labeled priority. Its scaling/zero-reference tests pass in serial, MPI and OpenMP; reverse derivatives match finite differences, and both new SA/SST MPI regressions pass. CI for the latest #2956 head is running. #2978 passes formatting; the regression workflow only triggers for PRs targeting develop/master, so its full matrix will run after retargeting it to develop following this merge.

ok, so this one goes first after updating regressions, then #2978

@bigfooted correction to my previous comment: the van Albada change in eef97f4 also changed the flow solver's VAN_ALBADA_EDGE, because the SIMD flow kernels use the same helper. 30P30N and dsma661 stopped converging, and updating their references in 48e871b hid that. This is fixed in 75621f9. The flow keeps its original limiter, and only the turbulence and species variables use vanAlbadaFunction with EPS. The V&V and axial_stage_2D references are back to their original values (30c52aa). The remaining species and rae2822_sst changes come from the boundary density upwinding. I will update those references from the next CI run.