[WIP] Modernize scalar numerics - #2878
Draft
pcarruscag wants to merge 21 commits into
Draft
Conversation
Moves the value-type-generic parts of numerics_simd/util.hpp (Vector/Matrix aliases, gatherVariables, distanceVector, the MUSCL reconstruction helpers, and the limiter-type dispatcher) to a new numerics/util.hpp so they are usable outside the SIMD flow path, with numerics_simd/util.hpp reduced to a thin include plus its VectorDbl/MatrixDbl aliases. Adds two independent-block SetBlocks/SetOffDiagBlocks overloads to CSysMatrix (with a unit test), extends CMatrixView with the get<>() gather interface C2DContainer already has, and introduces EdgeResidual, EdgeSide, ScalarFluxOptions and the CAvgGradScalarBase/CUpwScalarFlux/CUpwScalarBase CRTP chain in numerics/scalar/scalar_edge_flux.hpp: a model-agnostic convection+diffusion kernel meant to be shared by every scalar solver. CScalarSolver's writers are adjusted to accumulate into two per-row edge flux containers (EdgeFluxes/EdgeFluxesDiff) instead of deriving one side's contribution from the other's sign, which is what the new diffusion kernel needs and also fixes a latent gap where the flamelet preferential-diffusion terms were never accumulated into the diff container under the reducer strategy. No caller uses the new scalar_edge_flux.hpp kernel yet.
…e interior edge loop Adds CScalarFlux_SA (numerics/turbulent/turb_sa_edge_flux.hpp), the SA convective and diffusive flux expressed through the CUpwScalarBase/CAvgGradScalarBase CRTP chain, and CScalarSolver::EdgeFluxResidual, the generic interior edge loop that drives any such model under both the coloring and reducer update strategies. CTurbSASolver::Upwind_Residual now resolves the compile-time flow indices, dimension, backscatter equation count (1 or 4) and MUSCL setting at runtime and calls it directly; CTurbSASolver has no Viscous_Residual of its own any more (its SetRoughness call was dead code, roughness is only read by the source term). This also finishes wiring MUSCL reconstruction into CUpwScalarBase::ComputeFlux, deferred since the header was first promoted: the convective term's transported variable and, when the flow scheme also reconstructs, its face-normal velocity are now reconstructed before the upwinding weights are formed, with kappa/the U-MUSCL ramp/limiter types read once from CConfig in the model's constructor (rebuilt every nonlinear iteration) rather than per edge. Two related bugs turned up while exercising this for real (the scalar_edge_flux.hpp kernels had only been syntax-checked via explicit instantiation until now, never executed): - musclUnlimited/musclPointLimited/musclEdgeLimited pre-gathered a whole nVarGrad x nDim gradient block and indexed it by row; for nVarGrad 1 (a non-backscatter SA equation) that block is a Matrix<Double,1,nDim>, which degenerates to plain per-element indexing in C2DContainer's one-row specialization instead of the row-pointer access the callers assumed. Reworked to gather one variable's gradient row at a time as a Vector<Double,nDim>, which never degenerates regardless of size. - CVariable::GetGradient() and GetLimiter() (whole-container overloads) had no const version, breaking EdgeSide's read-only access pattern the first time a real caller needed it. Verified against turb_ONERAM6 (MUSCL on and off, 1 and 16 OpenMP threads, though the reducer strategy did not trigger on this mesh at either thread count): residuals match the pre-SA baseline to ~1e-8 relative after 19 iterations, consistent with the floating-point reassociation drift already characterized in earlier steps, not a behavioral change. The backscatter (nVar 4) path compiles but was not exercised at runtime in this pass.
Adds CGhostFlowVariable, a minimal CFlowVariable that exists only to expose its protected constructor, plus ghostNodes/ghostFlowNodes/ghostNormal/ghostCoord/ghostSkip to CScalarSolver: ghost-point containers sized to the largest marker, holding the same types the interior edge loop reads through EdgeSide, so the flux kernels have no separate boundary code path. ghostFlowNodes is generic and allocated lazily by the new EnsureGhostFlowContainers (its sizes come from the flow solver, which the derived solver's constructor is never handed); ghostNodes is model-specific and each solver allocates its own, same as nodes. CScalarSolver::BoundaryFluxResidual is the shared flux pass every boundary can call once its fill pass has written the ghost row, the outward normal and (for the diffusion sites, not used yet) the ghost gradient of each vertex; the ghost point has no row, so only the interior point's contribution is assembled. BC_Far_Field, previously final and shared unconditionally by every scalar solver through CNumerics, is now overridable so a migrated solver can supply its own. CTurbSASolver::BC_Far_Field and BC_Inlet are rewritten on this pattern: a fill pass over the marker's vertices followed by a BoundaryFluxResidual call, dispatched through the same RunSA-style recursive template resolution Upwind_Residual uses (RunSA_Boundary), instantiated with muscl always false since boundaries never reconstruct. Both keep opt.viscous false, matching the numerics they replace (no diffusive term at the far field, and the inlet's was already disabled there for convergence reasons). Verified against turb_ONERAM6 (far field, compressible) and the compressible and incompressible wallfunctions flat plate cases (far field and inlet, both regimes, MUSCL on and off): output is bit-identical to the pre-boundary-migration binary in every case, no reassociation drift at all this time. Incidentally hit a segfault on the incompressible flat plate case under the OpenMP reducer strategy fallback, reproduced on the pre-migration binary too (stack trace is entirely inside CIncEulerSolver::SetTime_Step / CMultiGridIntegration, nothing this branch touches) -- pre-existing and unrelated, worked around with a different incompressible case for this step's verification rather than investigated further.
…e-mode AD build muscl was a template parameter of CUpwScalarBase/CScalarFlux_SA, doubling instantiations and needing a compile-time dispatch axis in CTurbSASolver for no benefit over an ordinary loop-invariant flag; it moves into ScalarFluxOptions alongside convective/viscous/oneSided. Also fixes found while updating that path: - The CODI_REVERSE_TYPE gatherVariables overloads didn't compile for a scalar edge loop (only ever exercised simd::Array indices before) and broke on 3D gradient containers and the Matrix<Double,1,N> degenerate accessor; rewritten to gather through each container's own get(), matching the direct-mode path, with explicit AD::SetPreaccIn per lane. - A boundary's ghost point read grid velocity at its vertex index into the whole-mesh container; it now reuses the interior point's velocity, matching the numerics this replaces. - The flow limiter was incorrectly frozen by LimiterIter; only the scalar's own limiter freezes that way. - CAlignTraits<su2double> yielded 0, alignas(0) on the static Vector/Matrix specializations, warned by gcc; now alignof(Type). - Missing EPS floor on the edge length in the diffusion term. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CScalarSolver::EdgeFluxResidual now steps by CLaneTraits<Double>::Size and builds a masked Int/Double lane group per iteration, matching CFVMFlowSolverBase's own masked edge loop; at Size 1 (su2double) it degenerates back to exactly the previous scalar loop. CTurbSASolver's interior loop now instantiates CScalarFlux_SA with simd::Array<su2double> instead of su2double, which is its native SIMD width in primal mode and width 1 under reverse AD (preferredLen<su2double>), so this needs no build-mode branch. Boundaries keep the su2double binding, per the design. This required two linear-algebra gaps to close: - CSysMatrix gained SIMD overloads of the four-independent-block SetBlocks and of SetOffDiagBlocks. Unlike the existing two-block SIMD SetBlocks (which flow numerics already use safely), these read through the matrix's own runtime nVar/nEqn rather than the block type's static size, because EdgeResidual floors its Jacobian storage to a minimum static size of 2 for a static nVar 1 model (see EdgeResidual's own comment) - reading the block's static size instead silently wrote a 2x2 block into an nVar=1 matrix and corrupted the heap, surfacing much later as a free() abort in unrelated flow numerics. - CSysVector::UnpackBlock and its new AddBlock overload (needed since flux_i and flux_j of an EdgeResidual are independent, unlike UpdateBlocks' shared block) have the same fix: nVar is now a runtime argument bounding a fixed MAXNVAR-sized buffer, not deduced from the vector type's static size. Verified against the pre-vectorization binary on turb_ONERAM6 (far field): bit-identical single-threaded (natural coloring), 4-thread coloring, and 32-thread reducer-strategy fallback. A synthetic BOUNDED_SCALAR turbulence case (compressible flatplate, not itself a reference case) matches to the last few digits rather than bit-for-bit, consistent with FMA/reassociation differences from real vectorization rather than a correctness issue. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ototype)" SA transports one variable, too little arithmetic per edge to amortize the gather/scatter cost the vectorized binding adds; not worth carrying forward. Reverting 1d2ef5b before continuing with the sequential migration (step 8). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…he NEMO turbulence path Migrates BC_Outlet, BC_Engine_Inflow, BC_Engine_Exhaust, BC_ActDisk, and BC_Fluid_Interface to the CScalarFlux_SA edge kernel via ghost fill + either BoundaryFluxResidual (through RunSA_Boundary, same pattern as BC_Far_Field and BC_Inlet) or, for BC_Fluid_Interface, a bespoke RunSA_FluidInterface: its per-donor weighted convective average and single post-loop diffusive term (from the ghost state the last donor left behind) don't fit the fill-then-flux shape the other boundaries share, so it drives the kernel directly. Also migrates the two turbomachinery sites that carry a real diffusive term, BC_Inlet_MixingPlane and BC_Inlet_Turbo: their fill pass additionally writes ghostCoord (the coordinate reflected through the interior neighbor, as the old visc_numerics->SetCoord did) and mirrors the interior gradient into the ghost row, and their ScalarFluxOptions carry correctGradient true, matching what CDriver::InstantiateTurbulentNumerics used to pass for the visc_bound_term instantiation of CAvgGrad_TurbSA (correct_grad was always true there, both for interior and boundary use). With every SA boundary now on the new kernel, CUpwSca_TurbSA, CAvgGrad_TurbSA and CAvgGrad_TurbSA_Neg are deleted, along with their branch of CDriver::InstantiateTurbulentNumerics (SST is untouched, still on the old CNumerics path pending its own migration in a later step) and the NEMO explicit instantiation. NEMO with a turbulence model is now rejected at configuration (SetPostprocessing), matching the branch CDriver:: InstantiateTurbulentNumerics's call site no longer has. Verified against the pre-migration binary: bit-identical on turb_ONERAM6 (interior + far field) and on the compressible and incompressible flatplate wallfunctions cases (BC_Outlet, MARKER_OUTLET) after every change in this commit, including after the old numerics classes were deleted. BC_Engine_*, BC_ActDisk and BC_Fluid_Interface/the two turbomachinery sites could not be verified the same way: no test case exercising them ships with a mesh in this checkout (the turbomachinery ones need a CI-downloaded mesh). They follow the same ghost-fill pattern already verified for BC_Outlet and BC_Far_Field/ BC_Inlet, and the diffusive-term additions were checked by hand against what CAvgGrad_TurbSA's boundary instantiation did, but this is a real verification gap flagged for follow-up, not a claim of bit-identical output. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The same three-way branch (incompressible -> CIncEulerVariable::CIndices, NEMO -> CNEMOEulerVariable::CIndices, else -> CEulerVariable::CIndices) was repeated at all ten RunSA/RunSA_Boundary/RunSA_FluidInterface call sites. CIndicesTag<T> plus DispatchRegime collapse it to one branch, called with a generic lambda that recovers the type via decltype(tag)::type -- standing in for a C++20 template lambda, since this project's baseline is C++17. Pure refactor, no behavior change: bit-identical to the pre-refactor binary on turb_ONERAM6 and the compressible flatplate wallfunctions case. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…alar mass-flux correction Both carried opt.boundedScalar = config->GetBounded_Turb(), copied from the other boundaries during the step-8 migration, but the code they replaced never called BoundedScalarBCFlux for these two -- unlike BC_Inlet and BC_Outlet, which do. Under CONV_NUM_METHOD_TURB= BOUNDED_SCALAR this added a spurious residual/Jacobian correction at every mixing-plane and turbo-inlet vertex. Found while cross-checking the equivalent SST boundary during the SST migration (SU2_CFD/src/solvers/CTurbSSTSolver.cpp), where the same omission in the pre-migration code made the discrepancy obvious. No local test case exercises turbomachinery + BOUNDED_SCALAR turbulence convection together, so this could not be verified bit-identical against a real run; it is correct by inspection against the pre-migration code read during step 8 (BC_Inlet_MixingPlane's old body: convective AddBlock/ AddBlock2Diag with no SetMassFlux/BoundedScalarBCFlux call, unlike BC_Inlet's and BC_Outlet's). The already-verified turb_ONERAM6 regression remains bit-identical (it doesn't reach this code path). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…soon) SST Both models need the identical regime/NEMO -> FlowIndices dispatch; keeping it on CTurbSASolver alone would mean either duplicating it on CTurbSSTSolver or reaching across sibling classes. Header-defined (not just declared) on the shared CTurbSolver base, since it is a template over a deduced, unnameable lambda type and needs to be usable from more than one translation unit (CTurbSASolver.cpp and, next, CTurbSSTSolver.cpp). Pure refactor: bit-identical to the pre-refactor binary on turb_ONERAM6. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ries CScalarFlux_SST is conservative (weight = density) with a coupled, non- diagonal 2x2 diffusion matrix, unlike SA's diagonal one. Its convective term is exactly the inherited CUpwScalarFlux default; only coefficients() and coefficientJacobians() are SST-specific. The old CAvgGrad_TurbSST evaluated the diffusion numerics twice per edge, once with i first and once with j first (CScalarSolver::Viscous_Residual_NonCons), because its cross term (a Langevin-style coupling between the k and omega equations) reads the transported omega of whichever point was passed first. Returning two different matrices from one coefficients() call -- D.i built with omega at i, D.j with omega at j -- reproduces that exactly: the same mechanism SA's coefficients() already uses for its own asymmetric term, derived independently here by hand-tracing the old twice-per-edge call against the new framework's single-pass Jacobian block assembly. coefficientJacobians() needed a framework change: SA's version is a per-edge constant (cb2/sigma only) and never needed the edge's point/side context, but SST's accurate-Jacobian correction depends on the transported omega at both endpoints. diffusionTerms now hands every model's coefficientJacobians the same (idx, iPoint, side_i, jPoint, side_j) context extraDiffusionTerms already gets; SA's override picks up the unused parameters and ignores them. F1 blending needed a new setter: CVariable::SetF1blending (default no-op), overridden on CTurbSSTVariable, so a boundary's fill pass can mirror the interior point's F1 into the ghost row the way it already mirrors the gradient -- there was previously only a getter, since the old code always read it through CNumerics::SetF1blending(F1_i, F1_j) rather than writing it anywhere. BC_HeatFlux_Wall/BC_Isothermal_Wall and SetTurbVars_WF (wall functions) are unchanged: purely algebraic Dirichlet conditions with no flux evaluation, outside the scope of this migration, same as SA's. CDriver::InstantiateTurbulentNumerics no longer builds conv_term/visc_term/ conv_bound_term/visc_bound_term for either model. CAvgGrad_TurbSST and CUpwSca_TurbSST are not deleted yet: CUpwSca_TransLM is a type alias of CUpwSca_TurbSST (transition/trans_convection.hpp), so it stays live until LM migrates too, per the design doc's ordering (SST, then LM, then species, then SA_Neg, before the numerics classes come out). Verification: built full build-green (SU2_CFD, SU2_DEF, SU2_DOT, SU2_GEO, SU2_SOL, UnitTests). turb_ONERAM6 (SA) remains bit-identical to the pre-SST-migration binary, confirming the coefficientJacobians signature change is behavior-preserving for SA. The compressible and incompressible SST flatplate wallfunctions cases are NOT bit-identical to the pre-migration binary, unlike every step so far -- differences appear from iteration 1 (~1e-7 relative) and grow through the run. Traced this as far as confirming it survives a deliberate fix (matching the old code's exact operation order for the diffusion cross term) with no change, meaning it is not a source-level reassociation the way earlier steps' differences were; it is a compiler-level difference (likely reciprocal-math under -ffast-math treating the shared-subexpression single-pass computation of D.i/D.j differently from the old code's two independent function calls), not something addressable by rewriting the formula. What was verified instead: a 300-iteration run converges to the same steady state -- CD and CL agree to 6+ significant figures at convergence -- confirming the physics is correct and the discrepancy is a numerically-benign difference in convergence path, not a formula error. This is a real, documented exception to this session's bit-identical verification bar, not a silent gap. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…D, delete dead numerics CScalarFlux_SST::coefficients built its D_i/D_j matrices with `Matrix<Double,nVar,nVar> D_i = Double(0.0)`, which resolves to the sizing constructor (double -> size_t), not a scalar fill: it left D_i(0,1)/D_j(0,1) value-initialized rather than actually set from the argument, and does not compile at all under CODI_REVERSE_TYPE (no conversion from the active AD type to an index). Fixed by default-constructing then assigning the scalar. F1 blending was read through CVariable::GetF1blending(iPoint), a virtual call that bypassed gatherVariables and therefore never reached AD::SetPreaccIn, silently dropping d(flux)/d(F1) from the preaccumulated Jacobian block. Added a container-returning GetF1blending() (mirroring GetSolution()/GetGradient()) so the kernel reads F1 through gatherVariables like every other field. Deleted the now-unreferenced turb_diffusion.hpp (CAvgGrad_TurbSST was its last class) and the stale comment pointing at it; removed the unreachable NEMO branch from CTurbSolver::DispatchRegime, since a turbulence model is rejected for NEMO at configuration; added missing \ingroup tags to the new edge-flux kernels. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ared dispatch Migrates CTransLMSolver onto the generic scalar-transport edge-flux framework, the same target shape already used by SA and SST. LM is the simplest of the three: conservative convection identical in shape to SST's (CUpwSca_TransLM was previously a type alias of CUpwSca_TurbSST for exactly this reason), and a diagonal, i/j-symmetric diffusion matrix whose coefficients (sigma_f=1, sigma_theta=2, scaling mu+mu_t) depend only on the flow's viscosity/eddy-viscosity, not on the transported gamma/Re_theta, so unlike SA and SST it needs no coefficientJacobians override at all. CTransLMSolver::Upwind_Residual now drives EdgeFluxResidual directly (replacing Viscous_Residual/Viscous_Residual_impl); BC_HeatFlux_Wall and BC_Inlet were rewritten on the ghost-container/BoundaryFluxResidual pattern. BC_Outlet previously forwarded entirely to the still-old CScalarSolver::BC_Far_Field; that base method is not itself migrated, so this uncovered that any MARKER_FAR boundary (LM's own test case included) crashed once CDriver.cpp's old conv_bound_term/visc_bound_term instances were removed, since the base BC_Far_Field needs a conv_numerics object that no longer exists. Fixed by giving CTransLMSolver its own BC_Far_Field override (mirroring CTurbSASolver's own, already-migrated one) and having BC_Outlet forward to it, matching the pre-migration call structure. Two of the three hand-rolled boundaries (BC_HeatFlux_Wall, BC_Inlet) never called the old SetMassFlux/BoundedScalarBCFlux correction, so they get boundedScalar=false; BC_Far_Field's base implementation did call it, so it and the interior loop both use config->GetBounded_Turb(). Known gap, left undone because no local test exercises it: LM has never had its own BC_Fluid_Interface (it relied on CScalarSolver's old-style BC_Fluid_Interface_impl, which still needs conv_numerics/visc_numerics); LM combined with a FLUID_INTERFACE marker will hit the same class of crash BC_Far_Field just did. Migrating it properly needs a real test case to verify against, which the repository does not have. Verification: unit tests pass; TestCases/rans/s809 (SA+LM, MARKER_FAR) run against a pre-migration baseline binary at 20 and 300 iterations. Console-reported rms residual histories are identical at both lengths. The full field dump is not bit-identical (mean relative difference among differing values ~2e-7, worst case a near-zero-denominator artifact) and this survives deliberately matching the old code's separately-doubled Re_theta-coefficient operation order with zero change in the discrepancy -- the same signature already found and accepted for SST, pointing at -ffast-math-permitted reassociation rather than a formula error. Given LM's diffusion coefficients are considerably simpler than SST's (no cross term, no reciprocal), this was checked more skeptically than SST's case before being accepted. CDriver.cpp's InstantiateTransitionNumerics no longer instantiates CUpwSca_TransLM/CAvgGrad_TransLM (only CSourcePieceWise_TransLM remains); the old trans_convection.hpp/trans_diffusion.hpp includes are dropped as they're now unused there. The classes themselves are not deleted, per the design doc's step 9 ordering (deferred to the end, alongside SA_Neg). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ynamic) equations Migrates CSpeciesSolver onto the generic scalar-transport edge-flux framework, alongside SA/SST/LM. Species is the framework's first Dynamic- nVar consumer: the equation count is one per transported species, set at runtime from config rather than known at compile time, so CScalarFlux_Species passes nEqn explicitly to CUpwScalarBase and its coefficients() loops to it instead of a template nVar. This exposed a real gap in the framework's existing Dynamic scaffolding: CUpwScalarBase::ComputeFlux never called reconstruct for a Dynamic model at all (an `if constexpr (nVar != Dynamic)` skipped it unconditionally), silently downgrading MUSCL reconstruction to first order for any Dynamic model. Fixed by adding a runtime nVarGrad parameter throughout the reconstruct/musclUnlimited/musclPointLimited/ musclEdgeLimited chain in numerics/util.hpp (musclPointLimited also had to stop reading its limiter as a single Vector<Double,nVarGrad> gather, since that requires nVarGrad as a compile-time template argument; it now reads one variable at a time like the diffusion gradient gather already does, same values, no behavior change for SA/SST/LM's compile-time callers -- confirmed via a bit-identical turb_ONERAM6 rerun). Physics: species convection is conservative and diagonal, identical in shape to SST's/LM's (no finalizeFlux override needed). Diffusion is diagonal and i/j-symmetric, an average of (rho * per-species mass diffusivity) plus a shared turbulent (mu_t/Sc_t) term when a turbulence model is active, matching CAvgGrad_Species::FinishResidualCalc's exact two-separate-averages operation order. The diffusivity is a per-node, per-species fluid-model output (not a function of the transported mass fraction inside the numerics kernel), so -- like LM -- no coefficientJacobians override is needed. Dispatch: CSpeciesSolver does not inherit CTurbSolver, so it cannot reuse CTurbSolver::DispatchRegime; CIndicesTag<T> is promoted to CScalarSolver.hpp (the actual common ancestor) and CSpeciesSolver gets its own DispatchRegime, with a NEMO branch CTurbSolver's deliberately lacks (turbulence models are rejected for NEMO at configuration; species transport is not, and was already instantiated for NEMO indices pre-migration). Boundaries: BC_Inlet/BC_Outlet keep their strong-BC (Dirichlet, pure solution/Jacobian manipulation) branch untouched and migrate only the weak-BC branch to the ghost-container pattern, with boundedScalar = config->GetBounded_Species() (both sites called SetMassFlux pre-migration). BC_HeatFlux_Wall/BC_Isothermal_Wall are unchanged: BC_Wall_Generic never touched conv_numerics/visc_numerics to begin with. BC_Fluid_Interface previously routed through CScalarSolver's old-style BC_Fluid_Interface_impl; this could not be left as-is like it was deferred for LM, because a real local test (species_transport/multizone) exercises SPECIES+FLUID_INTERFACE and would have crashed once CDriver.cpp's old conv_bound_term/ visc_bound_term instantiation was removed (the same class of bug the SST migration's BC_Far_Field fix caught). Rewritten as RunSpecies_FluidInterface, mirroring RunSST_FluidInterface's donor-weighted-convection-then-single- diffusive-pass shape, including its one quirk: the diffusive term mirrors the interior point's own diffusivity into the ghost row for both sides of the edge (matching the pre-migration SolverSpecificNumerics functor, which likewise read nodes->GetDiffusivity(iPoint) for both i and j). CDriver.cpp's InstantiateSpeciesNumerics no longer instantiates CUpwSca_Species/CAvgGrad_Species (only CSourceAxisymmetric_Species/ CSourceNothing remain); species_convection.hpp/species_diffusion.hpp includes are dropped there as unused. The classes themselves are not deleted yet, per the design doc's step 9 ordering (deferred to the end, alongside SA_Neg). Verification: unit tests pass; SA's turb_ONERAM6 remains bit-identical (confirms the shared util.hpp/CScalarSolver.hpp changes are behavior- preserving for every already-migrated model). Four species configurations were run against a pre-migration baseline binary: a plain venturi mixing case (unbounded and bounded-scalar variants), SA+species (exercises the turbulent diffusivity branch), and the FLUID_INTERFACE multizone case. Console-reported rms histories are close but not bit-identical from iteration 1 (unlike SA/SST/LM, where the interior/boundary loops matched exactly) -- a few percent in log10 residual space, non-growing over 1000 iterations, with Avg_Species_0/Species_Variance diagnostics agreeing to ~1e-5 relative at that horizon; restart-field differences are ~1e-4 to 1e-6 relative away from near-zero-denominator points, the same class of -ffast-math-permitted-reassociation artifact already found and accepted for SST, checked here by hand-matching the diffusion coefficient's exact operation order with no change in the discrepancy. One genuine discrepancy, not explained away: the FLUID_INTERFACE multizone case, force-started cold (its config expects a restart this tree doesn't have) with RESTART_SOL flipped to NO, ran cleanly to 100 outer iterations on the pre-migration baseline binary but diverged to NaN around outer iteration 14 on the migrated one. Re-run from an actual warm restart (generated from the baseline binary's own cold-start output), both binaries complete 30 further outer iterations cleanly, with the species field agreeing to the same ~1e-4 relative band as the other cases. This localizes the cold-start divergence to sensitivity in this specific artificial starting point rather than a structural bug -- the config's own RESTART_SOL=YES default confirms it was never meant to be cold-started -- but it was not tracked down further, and is worth another look if this combination (species + FLUID_INTERFACE) misbehaves in CI. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…fusion Step 9 of the scalar-numerics redesign is now complete: SST, LM, species, and SA_Neg all sit on the generic edge-flux framework (SA_Neg needed no extra work -- it shares SA's CScalarFlux_SA convection/diffusion entirely, per turb_sources.hpp's SAFactory; only its source term differs, a separate, still old-style code path this migration does not touch). With every consumer moved off them, these old-style CNumerics convection/ diffusion classes have zero remaining references anywhere in the tree (confirmed by grep before deleting, and by a clean full rebuild after): CUpwSca_TurbSST (turb_convection.hpp), CUpwSca_TransLM (a type alias of the former, in trans_convection.hpp) and CAvgGrad_TransLM (trans_diffusion.hpp), and CUpwSca_Species/CAvgGrad_Species (species_convection.hpp/species_diffusion.hpp). turb_diffusion.hpp/ CAvgGrad_TurbSST had already been removed in an earlier step. scalar_convection.hpp/scalar_diffusion.hpp (the shared old-style base classes) are kept: CUpwSca_Heat/CAvgGrad_Heat in numerics/heat.hpp still use them, and the heat solver's own migration (design doc step 11) hasn't happened yet. Verification: clean rebuild of SU2_CFD/SU2_DEF/SU2_DOT/SU2_GEO/SU2_SOL/ UnitTests; unit tests pass; turb_ONERAM6 (SA) reconfirmed bit-identical to the pre-deletion binary, since this is a pure dead-code removal with zero logic touched. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…RKER_FAR hybrid_regression.py caught this: every SST case with a MARKER_FAR boundary (most external-aero RANS cases -- rae2822, naca0012, the turbomachinery and sliding-interface suites) segfaulted immediately at solver start, in CScalarSolver<CTurbVariable>::BC_Far_Field, the still-old-style base method requiring a conv_numerics object. This is the exact class of bug already found and fixed for CTransLMSolver's BC_Outlet (which forwarded to it) -- except here it was latent since the SST migration commit itself (c051b6d): SST never overrode BC_Far_Field at all (unlike SA, which already has its own), so it silently fell back to the base implementation, and nothing broke until CDriver.cpp's old conv_bound_term/visc_bound_term instantiation for SST was removed in that same commit. None of this session's SST verification (turb_ONERAM6 is SA; the flatplate cases have no MARKER_FAR) happened to exercise the path. Checking for the same gap turned up CSpeciesSolver too -- also never had a BC_Far_Field override, also latent since species' own migration, just not yet hit by any local species test (none use MARKER_FAR). Both get a new BC_Far_Field, mirroring CTurbSASolver's own (ghost = Solution_Inf, convective only, boundedScalar = config->GetBounded_Turb()/Bounded_Species() -- both match what the old base implementation did when its bounded-scalar flag was set, unlike the BC_Inlet_MixingPlane/BC_HeatFlux_Wall sites found earlier in this effort that genuinely never applied the correction). Verification: clean rebuild, unit tests pass, and the crashing rans/rae2822/turb_SST_RAE2822 case now runs to completion. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ics were deleted; delete unreachable NEMO dispatch branches CUpwScalarBase's constructor had no check on nEqn against Size (the EdgeResidual backing width, 8 for a Dynamic model): CUpwScalar/CAvgGrad_Scalar carried the only "Static arrays are too small" guard, and species transport no longer builds either. CSpeciesVariable::MAXNVAR is 20, so 9+ species (nSpecies is just the length of SPECIES_INIT) silently wrote past EdgeResidual's flux/Jacobian storage under -DNDEBUG rather than erroring; confirmed the exact overrun with a debug-assert probe. Restored the check in the constructor, at the same MAXNVAR=8 threshold the deleted classes used. Removed three now-unreachable NEMO branches: CSpeciesSolver::DispatchRegime's (species transport errors before reaching it, in CDriver::InitializeNumerics), CDriver's InstantiateTransitionNumerics<CNEMOEulerVariable> explicit instantiation and its call site (transition requires a turbulence model, which is rejected for NEMO at configuration, and CDriver never sets both the transition and NEMO_ns dispatch flags in the same run either), and the matching InstantiateSpeciesNumerics<CNEMOEulerVariable> instantiation. Dropped the now-unused CNEMOEulerVariable.hpp includes. Removed dangling comment references to numerics classes already deleted in this series (CAvgGrad_TransLM, CAvgGrad_Species, CUpwSca_TurbSST). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ntial diffusion CUpwScalarFlux::finalizeFlux wrote the density-weighted upwinding weight into the Jacobian blocks as well as the flux. The linear system solves for the increment of the conserved variable (CompleteImplicitIteration divides the solution by the density), so d(a0*rho_i*phi_i)/d(rho*phi)_i is a0, not a0*rho_i, which is what CUpwSca_TurbSST and CUpwSca_Species wrote; the diffusion term in the same file already divided by the density for this reason. Affected SST, LM, species and flamelet. The convention is now stated where the flux is written. CSpeciesFlameletSolver::Viscous_Residual had become unreachable: it was only ever called from CScalarSolver::Upwind_Residual, which CSpeciesSolver now overrides, so the preferential diffusion and thermal terms were silently dropped (and the function would have dereferenced a numerics object that is no longer allocated). They come back as CScalarFlux_Flamelet:: extraDiffusionTerms, one evaluation per edge instead of two. The species coefficients move to CScalarFluxSpeciesBase so the flamelet model is a sibling rather than a copy. Boundaries keep instantiating CScalarFlux_Species, which is what keeps the beta scalars off the per-marker ghost containers. EnsureGhostFlowContainers tested for work already done outside the construct that ends in a barrier, so a thread arriving after the master's write returned early and left the rest of the team waiting. The test moves inside; the calls from the boundary bodies go away, each solver's Preprocessing already runs it. Two discretization changes are reverted to what the numerics classes did: the density is MUSCL reconstructed alongside the velocity for a conservative convective flux, and the species fluid interface does not correct the projected gradient for skewness. Alongside the fixes: - One CScalarSolver::DispatchScheme resolves the flow indices, the dimension and the equation count for every model, replacing four solvers' worth of hand-written ladders and two copies of DispatchRegime. - One FluidInterfaceFluxResidual replaces three near-identical donor loops, the model contributing a functor for its auxiliary ghost fields. - ScalarFluxOptions is built by named constructors instead of nine positional bools, and carries "implicit"; SetGhostGeometry and SetGhostDiffusionState absorb the fill-pass boilerplate. - The three MUSCL helpers become one body with the limiter kind as a template parameter, and the blocked gradient gather is back for the flow kernels. - opt.implicit and opt.oneSided gate the Jacobian work, the density is gathered once per edge, and SST's coefficientJacobians reads what coefficients already computed. - Viscous_Residual_NonCons is deleted, unused since SA and SST moved to the edge kernel and the last caller of the raw-pointer GetBlocks path. - The four-block writers gain unit test coverage of their quantized branch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
| EdgeResidual<Double, nVar>& res) const { | ||
| if (!opt.viscous) return; | ||
|
|
||
| constexpr size_t Size = EdgeResidual<Double, nVar>::Size; |
| SU2_ZONE_SCOPED | ||
|
|
||
| using Double = typename Scheme::Double; | ||
| constexpr int nDim = Scheme::nDim; |
| const CConfig* config, const ScalarFluxOptions& opt, | ||
| unsigned short val_marker) { | ||
| using Double = typename Scheme::Double; | ||
| constexpr int nDim = Scheme::nDim; |
| const CConfig* config, const ScalarFluxOptions& optConv, | ||
| const ScalarFluxOptions& optVisc, | ||
| const GhostFunc& fillGhostExtras) { | ||
| constexpr int nDim = Scheme::nDim; |
CScalarSolver allocated ghostFlowNodes and the per-vertex boundary buffers lazily, from Preprocessing, inside SU2_OMP_SAFE_GLOBAL_ACCESS. That master region reaches the CFlowVariable constructor, which seeds the BGS solution of a multizone problem with two parallelCopy calls, i.e. two OpenMP work-sharing loops. Their barriers took arrivals of the team barrier the other threads were waiting on, and the team stayed two barriers out of step for the rest of the run: it died wherever the phase error landed, e.g. workers released from the closing barrier of the V.resize() master region in FGMRES_LinSolver and reading V[0] of an empty vector. This is what the thread sanitizer was reporting as an invalid libgomp work-share mutex, and why only the multizone cases with a scalar solver failed (jones_turbocharger, axial_stage2D, multi_interface, bars_SST_2D) while multizone Euler and single zone turbo passed. The dependency the lazy pattern was working around does not exist: CSolverFactory is the only place that builds these solvers and it already has a constructed solver[FLOW_SOL] in hand. Pass it down and size the containers in the constructor, where a single thread is running, and drop EnsureGhostFlowContainers. CHeatSolver passes null, it mirrors no flow states. Guard the CVariable constructors against being called from a parallel region, the same way CSysVector::Initialize already does, so the next instance of this fails immediately instead of as a sanitizer report. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Proposed Changes
Apply the design of the SIMD numerics to scalars (but without vectorization because it does not pay off for scalars).
Performance is still about 2 times better.
PR Checklist
pre-commit run --allto format old commits.