pressure based corrections - #2896
Merged
Merged
Conversation
TestCases/incomp_navierstokes/lid_driven_cavity/incomp_liddrivencavity.cfg and incomp_pb_liddrivencavity.cfg were added but never referenced by any regression list, and no mesh for them (square_50x50.su2) exists anywhere in this repo or in the su2code/TestCases mesh repo. tutorials.py already exercises the same case (Lid Driven Cavity Flow, pressure-based) via ../Tutorials/incompressible_flow/Inc_Lid_Driven_Cavity/ incomp_pb_liddrivencavity.cfg, which lives in the external su2code/Tutorials repo referenced by this PR's description (Tutorials#86) and carries its own mesh. These two files were unrunnable duplicates of that. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
RMS_MASSFLUX is not a registered history output field in CFlowIncOutput, so it was silently dropped from SCREEN_OUTPUT rather than printed. This affected incomp_pb_NACA0012.cfg and incomp_pb_cylinder.cfg (the flat plate and lid-driven-cavity PB configs did not have this line). Replaced with RMS_PRESSURE, which for the pressure-based solver is already remapped in CFlowIncOutput::LoadHistoryData to the Poisson solver's residual, so this now shows real, useful convergence data instead of a silently-ignored field. The regression test_vals in serial_regression.py and parallel_regression.py are unaffected: TestCase.py matches the trailing N columns of the output row, and RMS_PRESSURE is inserted before the existing RMS_VELOCITY-X/-Y, LIFT, DRAG columns, so their relative order and count are unchanged. Verified by direct run against the existing registered values. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Space_Integration runs inside an active OpenMP parallel region (CSingleGridIntegration.cpp), so every thread enters Source_Residual independently. Two separate bugs meant its boundary and source-term contributions were wrong under threading: 1. The per-marker vertex loops (INLET_FLOW, FAR_FIELD, and the PRESSURE_OUTLET case of OUTLET_FLOW) had no #pragma omp for at all, so every thread executed the full vertex range redundantly and LinSysRes.AddBlock() applied each boundary mass-flux correction nThreads times over. Fixed by adding SU2_OMP_FOR_DYN(OMP_MIN_SIZE) / END_SU2_OMP_FOR around each loop, matching the existing CIncEulerSolver::BC_Far_Field / BC_Outlet idiom. 2. The interior edge loop wrote (edgeMassFluxes[iEdge] + MeanHbyA) directly into LinSysRes via AddBlock/SubtractBlock at both edge endpoints. Under the ReducerStrategy fallback (used when edge coloring efficiency is poor), edges sharing a point can be processed concurrently by different threads, racing on that point's LinSysRes row. Fixed by staging the combined per-edge value into a new EdgeSourceFlux member (write-once per edge index, race-free regardless of coloring) and scattering it into LinSysRes in a second, point-partitioned pass. Could not reuse the existing SumEdgeFluxes()/EdgeFluxes machinery here because it resets LinSysRes before scattering, which would erase the diffusion residual Viscous_Residual already assembled earlier in the same pass; the new scatter loop is additive instead. Verified against pb_rough_flatplate_incomp.cfg: on the unfixed code, at the default OpenMP thread count on a 32-core machine, the case diverges (Residual > 10^20) by inner iteration ~7. With only fix (1) applied it still diverges at OMP_NUM_THREADS=2, confirming (2) was independently necessary. With both applied, mpirun -n 2 + OMP_NUM_THREADS=2 and the hybrid_regression.py invocation (SU2_CFD -t 2) both reproduce the registered test_iter=10 values exactly. Also re-verified the existing hydrofoil and cylinder PB regression cases still match serial_regression.py/parallel_regression.py bit-for-bit with these fixes and the current develop merged in. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
hybrid_regression.py (SU2_CFD -t 2, single MPI rank) had zero PB coverage, despite being the only regression driver that exercises the Poisson solver's OpenMP path (the other regression scripts build with OpenMP disabled). Added the hydrofoil and flat-plate PB cases; picked these over the lid-driven cavity because every marker there carries zero boundary mass flux, which would make the now-fixed Source_Residual threading bugs invisible (the redundant nThreads-times accumulation of a zero correction is still zero). test_vals measured directly against SU2_CFD -t 2 with the OpenMP fixes applied; closely tracks (to 5-6 significant figures) the values already registered in serial_regression.py/parallel_regression.py for the same configs, as expected. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…phere (3D) cases Both are copies of existing density-based configs with KIND_INCOMP_SYSTEM= PRESSURE_BASED and the pressure-based solver options added, reusing the existing mesh in place. bend: the only CGNS mesh and the only strongly non-orthogonal internal geometry in the PB test suite. At the DB sibling's CFL_NUMBER=10 the case does not diverge outright within 10 iterations but does not converge either (rms[U] grows from -17.5 to +1.8, CD swings to -6654). Dropping to CFL_NUMBER=1.0 gives a clean, monotonically-converging trajectory. This 10x CFL sensitivity on the one non-orthogonal mesh in the suite is worth keeping in mind alongside the missing skewness correction in the pressure-correction Laplacian. sphere: the only 3D case in the PB test suite. Runs cleanly at CFL=10 with no CFL sensitivity issue. RESTART_SOL=NO unlike the DB sibling, since the pressure-based solver does not restore edge mass fluxes from a restart file. Registered in parallel_regression.py only, matching where the DB sibling lives. test_vals measured directly against the built solver (serial and np=2 for bend; np=2 for sphere, matching where each is registered). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
inc_flatplate_pb (pb_rough_flatplate_incomp.cfg) diverges to Residual > 10^20 around inner iteration 85-93 even under fully correct, non-threaded execution (verified against the -Dwith-omp=false build). The registered test_iter=10 checkpoint itself is unaffected - TestCase.run_test only greps the log for that one row and never checks the process exit code - but shipping coverage of a case known to diverge later is misleading as a regression fixture and not worth it until the underlying formulation issue (likely F1: the A_p artefact at strong-BC nodes, worst on this case's high-aspect-ratio near-wall cells) is addressed. inc_euler_naca0012_pb remains in hybrid_regression.py and already exercises the same fixed code paths (nonzero-flux INLET_FLOW and OUTLET_FLOW markers) without this problem - it runs to completion cleanly with no divergence at any point. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
SetMomCoeff, ComputeEdgeMassFluxesRhieChow, ComputeHbyA and ApplyPressureVelocityCorrection are called directly from CPBFluidIteration::Iterate, outside any parallel region. Their internal SU2_OMP_FOR_* pragmas were consequently orphaned no-ops (a work-sharing construct outside a parallel region just runs on the calling thread), so the entire pressure-correction half of the algorithm was single-threaded regardless of the requested OpenMP thread count - the "Poisson solver is a major bottleneck" behaviour flagged in the PR description. Wrapped each of the four call sites in its own SU2_OMP_PARALLEL block, matching the existing AdaptCFLNumber idiom later in the same function. Deliberately did not wrap the MultiGrid_Iteration/SingleGrid_Iteration calls that sit between them - those already open their own parallel region internally, and nesting one parallel region inside another would (with nested parallelism off, the default) collapse the inner region to a team of size 1, causing every outer thread to redundantly re-run the entire momentum or Poisson solve rather than sharing it. Wrapping the call sites exposed one previously-dormant, genuine correctness bug and several previously-dormant performance gaps, none visible before because the code only ever ran on the calling thread: - End of ApplyPressureVelocityCorrection: the loop accumulating the mass-flux correction back into EdgeMassFluxes had no work-sharing construct at all. Once the call site is genuinely parallel this is an accumulation executed redundantly by every thread for every edge, silently multiplying the correction by nThreads - the same class of bug as the Source_Residual fix. Added SU2_OMP_FOR_STAT. - The edge loops computing the Rhie-Chow mass fluxes and the mass-flux correction, and the OUTLET_FLOW/FAR_FIELD marker-vertex loops in ApplyPressureVelocityCorrection, had no work-sharing construct either. These are all safe-but-wasteful idempotent per-slot overwrites (not accumulations), so they were not correctness bugs, but every thread was doing the full amount of work redundantly. Added the appropriate SU2_OMP_FOR_STAT/SU2_OMP_FOR_DYN to each so they actually partition. - SetMomCoeff and ComputeHbyA already had correct SU2_OMP_FOR_STAT on their main loops; no changes needed there. Verified: the flat-plate, cylinder and hydrofoil PB regression cases still match their registered test_iter values after this change (both at -t 2 and at mpirun -n 2 with the default OpenMP thread count) - no correctness regression from either the new parallelism or the fixes it exposed. A 100-iteration timing run of the cylinder case went from 28.1s at -t 1 to 6.4s at -t 8 (about 4.4x), with the -t 1 result matching serial_regression.py's registered value bit-for-bit and -t 8 differing only at the 4th-5th decimal (expected floating-point summation-order noise from a different thread partition). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…er plumbing - CSolverFactory::CreateSubSolver: metaData.integrationType was set to SINGLEGRID unconditionally, even when the pressure-based solver is off and no CPoissonSolver is created. In practice this was harmless (CreateIntegrationContainer only queries metadata for non-null solvers, and the allocatedSolvers map is only populated when genericSolver != nullptr), but it read as a real bug and would break the moment either of those guards moved. Moved the assignment inside the same if as the allocation. - CScalarSolver::LoadRestart was changed from a pure virtual to an empty override, silently dropping the compile-time guarantee that every scalar solver implements restart loading, for every current and future subclass, to avoid four lines in CPoissonSolver. Restored `= 0` and added an explicit empty override in CPoissonSolver with a comment explaining why: the pressure correction is reset to zero every iteration (see Preprocessing), so it carries no state a restart needs to provide. Checked all four other CScalarSolver subclasses (CHeatSolver, CTurbSolver, CSpeciesSolver, and CTransLMSolver via CTurbSolver) already implement it, so this is a pure compile-time check with no other source changes needed. - CSolver::GetCommCountAndType: MPI_QUANTITIES::MOM_COEFF declared COUNT_PER_POINT = nDim for a quantity that packs and unpacks a single scalar (confirmed against the Initiate/CompleteComms pack sites, which only touch bufDSend[buf_offset]). Fixed to 1, removing a 2-3x oversized comm buffer for this quantity. Verified: incomp_pb_cylinder.cfg (mpirun -n 2) still matches parallel_regression.py's registered test_iter=10 values exactly after all three fixes. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…runs CIncEulerSolver::Source_Residual reads the pressure gradient (GetGradient_Primitive at prim_idx.Pressure()) to build the pressure-based solver's momentum source term (V * grad(p)). For a viscous run CIncNSSolver::Preprocessing computes this gradient unconditionally, but CIncEulerSolver::Preprocessing (the inviscid path, used directly by SOLVER=INC_EULER) only ever wrote into Gradient_Reconstruction, and only inside the MUSCL block. That happens to alias into Gradient_Primitive when the reconstruction and base gradient methods match (CFlowVariable's Gradient_Reconstruction is a reference to Gradient_Primitive unless GetReconstructionGradientRequired() is true), which is why incomp_pb_NACA0012.cfg has always worked - it sets MUSCL_FLOW=YES with both methods defaulting to GREEN_GAUSS. With MUSCL_FLOW=NO, or mismatched gradient methods, the gradient was never computed at all and the momentum equation silently lost its pressure forcing term. Added an unconditional gradient computation, gated on pressure_based, mirroring CIncNSSolver::Preprocessing's unconditional block (same GREEN_GAUSS/WEIGHTED_LEAST_SQUARES-only coverage - LEAST_SQUARES has the same gap here as everywhere else this pattern appears, tracked separately). Verified the fix actually matters, not just that it compiles: built and ran incomp_pb_NACA0012.cfg with MUSCL_FLOW=NO forced. Before this fix, CL swings from 1.2 to -28.9 and CD from 0.35 to -13.6 within 13 iterations - the SIMPLE algorithm becomes unstable almost immediately without the pressure forcing term. After the fix, the same case converges smoothly and monotonically. The registered incomp_pb_NACA0012.cfg regression case (which already had MUSCL_FLOW=YES) is unaffected - this now recomputes the same gradient the MUSCL block already computed via aliasing, bit-for-bit identical result at test_iter=20 - and the cylinder, bend and sphere PB cases still match their registered values exactly. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…d code) A prior review pass flagged CIncEulerSolver::ComputeEdgeMassFluxesRhieChow and CPoissonSolver::Postprocessing for handling GREEN_GAUSS and WEIGHTED_LEAST_SQUARES but not LEAST_SQUARES, reasoning that under LEAST_SQUARES the Poisson operator would assemble against a stale gradient (since CommonPreprocessing has the same gap) while the flux correction used the compact form, breaking mass conservation. That reasoning turns out not to matter: NUM_METHOD_GRAD=LEAST_SQUARES (unweighted, as opposed to WEIGHTED_LEAST_SQUARES) is already rejected unconditionally for every SU2 case by CConfig::SetPostprocessing (Common/src/CConfig.cpp, "LEAST_SQUARES gradient method not allowed for viscous / source terms") - confirmed by actually trying it against both a viscous and an inviscid pressure-based case; both abort at config validation before the solver ever runs. This restriction predates this PR and applies regardless of whether the pressure-based solver is in use, so the code path this review comment was worried about is unreachable dead code, not a live conservation bug. Added the LEAST_SQUARES case anyway, alongside my own newly-added gradient computation in CIncEulerSolver::Preprocessing (which has the identical characteristic) - harmless, and consistent with SetPrimitive_Gradient_LS/SetSolution_Gradient_LS themselves already being written to handle both weighted and unweighted internally. Verified no change in behavior: incomp_pb_cylinder.cfg and incomp_pb_NACA0012.cfg (both GREEN_GAUSS) still match their registered test_iter values bit-for-bit. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Mechanical, non-behavior-changing cleanup across the PB solver files, plus two correctness-adjacent fixes: - Ran the missing whitespace/EOF-newline cleanup pre-commit would have caught: stripped trailing whitespace and fixed missing final newlines in pressure_based.hpp, CPoissonSolver.cpp and CPoissonVariable.cpp. - Deleted CPBConvection_Base::MeanPressure (computed, never read). - Deleted the GetExtraOutput() block copied verbatim from CHeatSolver into CPoissonSolver's constructor - OutputHeadingNames was allocated but never populated, and nothing in CPoissonSolver ever reads OutputVariables/nOutputVariables. - Removed the dead PRef_Point/Pref_local lookup in ApplyPressureVelocityCorrection: both branches unconditionally left PCorr_Ref at 0.0, so the lookup had no effect. Replaced with a comment stating plainly that the reference pressure is fixed at 0 and explaining when that stops being sufficient (a fully closed domain has no pressure boundary condition and needs a real datum instead). - CDriver::InitializeNumerics declared pressure_based and poisson with identical initializers; kept pressure_based, updated the one use site. - CFluidIteration::CommonAuxiliarySolvers took output, surface_movement, grid_movement and FFDBox but used none of them; trimmed the signature and both call sites. - Fixed the \breif typo in SUB_SOLVER_TYPE::POISSON's comment. - Fixed misleading attribution and documentation left over from copy-pasting existing files as a starting point: CPBFluidIteration.*, CPoissonSolver.* and CPoissonVariable.* were credited to "F. Palacios, T. Economon" / "O. Burghardt" (the heat/turbulence solver authors whose files these were copied from) instead of this PR's actual author; CPBFluidIteration.hpp's class doxygen said "\class CFluidIteration"; CPoissonVariable.hpp's brief said "heat equation solver". - config_template.cfg documents default option values, but the PB block set KIND_INCOMP_SYSTEM=PRESSURE_BASED, USE_AUTOMATIC_RELAXATION_FACTORS=YES, PISO_CORRECTIONS=2 and POISSON_LINEAR_SOLVER_ITER=1000, none of which match the actual registered defaults (DENSITY_BASED, NO, 1 and 10 respectively per CConfig::SetConfig_Options) - anyone using the template as a starting config would silently get the pressure-based solver and different settings than the defaults it claims to document. Fixed all four, and PISO_CORRECTIONS's "= " spacing to match the file's KEY= value convention. (Checked the "missing % separator before SOLID ZONE HEAT VARIABLES" item from the same review pass - the section transition already matches the file's established pattern exactly; that finding was a mistake, left alone.) - PrepareImplicitIteration re-queried config->GetKind_Incomp_System() instead of using the cached pressure_based member; switched to the member. - BC_Inlet used the raw offset V_inlet+1 for velocity; changed to V_inlet+prim_idx.Velocity() for consistency with the rest of the function (prim_idx.Velocity() is 1 for this variable layout, so this is behavior-preserving). - Added a comment to CPoissonSolver::Source_Residual's boundary correction loop explaining why it negates via -= instead of negating Normal first (CIncEulerSolver's convention) - both are correct, the inconsistency was just undocumented. - CIncEulerVariable::strongBC was allocated (uninitialized) for every incompressible run, not just pressure-based ones. Guarded the allocation on KIND_INCOMP_SYSTEM==PRESSURE_BASED and initialized to false; verified every SetStrongBC/GetStrongBC/ResetStrongBC call site is already reachable only when pressure_based is true. One review finding from this same pass turned out to be wrong on closer inspection, corrected rather than silently dropped: CFlowIncOutput's two convFields.empty() checks are mutually exclusive on pressure_based (one checks it, the other checks its negation), so the second one does fire correctly when needed - it is not dead code. Rewrote as an if/else for clarity, which was the underlying suggestion worth keeping regardless of that mistaken premise. Verified no behavior change: all five pressure-based regression cases (hydrofoil, cylinder, flat plate, bend, sphere) still match their registered test_iter values exactly. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…essure-based solver None of these five are wrong for a particular mesh or marker set - they are wrong for every pressure-based run regardless of configuration, so they can be caught once at config-validation time instead of being audited case by case: - MGLEVEL > 0: the Poisson solver only ever runs on the finest grid (CSolverFactory already sets INTEGRATION_TYPE::SINGLEGRID for it); a coarse-level POISSON_SOL would be allocated numerics with nothing driving them. - TIME_DOMAIN= YES: the Rhie-Chow mass flux interpolation has no pseudo-transient term. - MATH_PROBLEM= (CONTINUOUS|DISCRETE)_ADJOINT: no adjoint of the segregated SIMPLE/PISO correction exists. - MARKER_PERIODIC: the marker switches in ApplyPressureVelocityCorrection and Source_Residual have no PERIODIC_BOUNDARY case. Previously this fell through to the generic "boundary condition not implemented" error, but only once the solver reached its first iteration; guarding it here fails at config time instead. - KIND_STREAMWISE_PERIODIC != NONE: SetBeta_Parameter is disabled under the pressure-based path, so the option would silently do nothing. Placed in CConfig::SetPostprocessing's existing "Check for unsupported features" section, after the two places that can still zero nMGLevels for other reasons (adjoint / TIME_STEPPING deactivation), so the MGLEVEL check only fires on multigrid that would actually reach the solver. Verified each guard fires with the intended message by flipping one option at a time on incomp_pb_NACA0012.cfg, and that all four currently-registered pressure-based regression cases (hydrofoil, cylinder, bend, sphere) still match their registered test_vals with zero delta under mpirun -n 2, and that an unmodified density-based incompressible case is unaffected. This does not attempt the rest of the compatibility audit (ActDisk, Riemann/Giles, energy equation, species, dynamic grid, etc.) - those need a per-case verdict, not a blanket rule, and are tracked separately. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…tum coefficient SetMomCoeff reads A_p = Jacobian(iPoint,iPoint)(1,1)/rho directly off the momentum Jacobian diagonal. At any point where a strong velocity BC applied Jacobian.DeleteValsRowi (walls, inlets, far-field inflow), that call has already zeroed the entire row and written 1.0 on the diagonal before PrepareImplicitIteration adds Vol/dt to every diagonal - so A_p there is (1 + Vol/dt)/rho, a row-deletion artefact with the wrong units, not a momentum coefficient. This isn't just a convergence-rate cost: A_p feeds Coeff_Mom = 0.5*(d_i + d_j) in Rhie-Chow on every edge touching a wall or inlet, so it changes the converged solution itself. SetMomCoeff now runs two passes instead of one. The first reads the raw A_p as before and flags a point whose momentum row has zero coupling to every neighbour - the exact signature DeleteValsRowi leaves, since a genuine interior or weak-BC row always has nonzero convective/diffusive neighbour coupling (this is an exact floating-point equality check against a value the code itself sets to precisely 0.0, not a tolerance heuristic). The second pass substitutes, at a flagged point, the average raw A_p of its non-flagged, locally-owned neighbours - extrapolating from the interior momentum operator rather than reconstructing it from scratch, since the operator itself only exists as an assembled sparse Jacobian with no closed form to re-evaluate outside the pass that already produced the (now BC-overwritten) matrix. Off-rank neighbours are skipped; if every neighbour is also flagged, the artefact value is kept rather than left undefined. Validated against the DB↔PB agreement criterion this card sets out to satisfy, not just a green test suite: ran incomp_pb_cylinder.cfg and its density-based sibling to actual Cauchy convergence (not test_iter) in a scratch harness. DB converges to CD=2.865931, CL=0.000068. PB before this fix converges to CD=2.800851 (2.27% off DB), CL=-0.001283. PB after converges to CD=2.852209 (0.48% off DB - a 4.7x reduction in disagreement), CL=-0.000974, closer to the physically correct near-zero value for symmetric flow past a cylinder. This does not fix the flat plate's divergence (still blows up around inner iteration ~95-105 with the fix applied, versus ~85-93 without it - a small delay, not a resolution), so A_p alone does not explain that case's instability; F3's CFL ceiling and F6's missing skewness correction remain open candidates given the mesh's face-area aspect ratio up to ~42858. Re-baselined test_vals for all four currently-registered pressure-based cases, since this genuinely changes the converged answer rather than just moving it around: inc_lam_cylinder_pb and inc_lam_bend_pb move measurably in both serial_regression.py and parallel_regression.py, inc_lam_sphere_pb moves moderately in parallel_regression.py, and inc_euler_naca0012_pb moves by ~1e-6 (its only strong velocity BC, MARKER_INLET, sits far upstream of the monitored airfoil surface). The bend case's large test_iter=10 delta was checked against a 100-iteration run and is a stable, monotonically converging trajectory, not a new instability - it settles at CD~=0.4525 versus the previous run's still-transient CD=1.88 at iteration 10. hybrid_regression.py's copy of the naca0012 case is unaffected (bit-for-bit identical) since it isn't re-registered here. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
ApplyPressureVelocityCorrection computed the automatic-relaxation coefficient as factor = sum over iDim in [0,nDim) of Jacobian.GetBlockView(iPoint,iPoint)(iDim,iDim), then alpha_p = (Vol/dt) / (factor + Vol/dt). Three things were wrong with it: - Off by one: block row/column 0 is the continuity (pressure) row, not a velocity direction. Starting the sum at iDim=0 mixed the pressure row's diagonal (set to exactly 1.0 by DeleteValsRowi at every strong-BC node) into what is supposed to be a momentum coefficient. - It summed every velocity direction's diagonal instead of taking one, so alpha_p came out dimension-dependent (2D and 3D would disagree on the same physics) - and in 3D it summed pressure + u-momentum + v-momentum while omitting w-momentum entirely (three terms, iDim in [0,3), off a 4-wide block for a 3D problem with the pressure row included). - Missing /rho: a_P = dR/d(rho u) has mass/time units, matching the Vol/dt term it's compared against requires the same /density SetMomCoeff already applies for exactly this reason. Replaced the loop with a single read of the x-momentum row, matching SetMomCoeff's own convention that this coefficient is the same in every direction: factor = view(1,1) / nodes->GetDensity(iPoint). USE_AUTOMATIC_RELAXATION_FACTORS is NO by default (config_template.cfg, fixed earlier this pass) and none of the four currently-registered pressure-based cases set it, so this branch was dead code in the existing regression suite - confirmed all four are bit-for-bit unchanged. To actually exercise the fix, added pb_sphere_urf.cfg (the sphere case with USE_AUTOMATIC_RELAXATION_FACTORS=YES; 3D because the bug is invisible in 2D, where the missing w-momentum term and the extra pressure-row term partially offset by coincidence of dimension) and registered it in parallel_regression.py. Measured the fix's actual effect with a temporary diagnostic print (verified, then removed - not part of this diff): before the fix, factor~=9.65, alpha_p~=0.161 at one interior point; after, factor~=4.38, alpha_p~=0.297 - nearly double, in the direction the bug analysis predicts, even though this case's density is 1.0 (so the missing /rho term alone has no numerical effect here - the off-by-one and summed-directions bugs are what this particular case exercises). The case now runs stably for its full 10 iterations with the fix. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…cobian MassFlux is set once per outer iteration from Rhie-Chow and held fixed during the momentum solve, so its exact derivative with respect to velocity is zero at this stage. ComputeJacobian's momentum block instead wrote J[i+1][j+1] = scale * rho * (v_velocity[i]*Normal[j] + proj_vel*delta_ij) The proj_vel*delta_ij term is the real d(m_f*v_adv)/du at fixed m_f and correctly carries the advection scheme's weight (upwind: 1 or 0; central: 1/2), since v_adv is what the scheme actually upwinds or averages. The other term instead treats the frozen mass flux as if it depended on the neighbouring velocities, weighted by the same upwind/downwind selector - which for CPBConvection_Upwind is 1 or 0, not the 1/2 a term standing in for something that (if it depended on velocity at all) would depend on it symmetrically. Tried reweighting that term to a fixed 1/2 first, as directly suggested by the analysis. It broke incomp_pb_NACA0012.cfg (the only registered case with MUSCL_FLOW=YES) catastrophically - CL/CD diverge past 10^17 by inner iteration 8 at the case's registered CFL_NUMBER=10, though it converges cleanly at CFL_NUMBER=1, confirming this was a stability-margin regression rather than a broken derivation, and the opposite of what a diagonal-dominance argument for this term would predict. Fell back to the alternative already available in the same analysis - drop the term entirely, since it has no legitimate derivative to represent once the mass flux is frozen - and verified: stable at CFL=10 for the hydrofoil, and all four other registered pressure-based cases (cylinder, bend, sphere, sphere+URF) run cleanly with plausible, different, non-crashing values. The most consequential effect: pb_rough_flatplate_incomp.cfg, which has diverged catastrophically since before any of this session's fixes (by inner iteration ~85-93 with nothing fixed, ~95-105 with F1 alone), now survives to roughly iteration 255 with F1+F3 both applied before eventually diverging - a ~2.8x extension, not a resolution. Not claiming this case fixed; it still fails eventually, just much later, and its already-registered test_vals at test_iter=10 are unaffected to 5-6 significant figures (this term's effect only shows up well into the case's trajectory). F6 (missing skewness correction) and F9 (hard 0.5 interpolation weights) are the leading remaining candidates given this mesh's face-area aspect ratio up to ~42858. Re-baselined test_vals for the four already-registered pressure-based cases that do move measurably at their checked iteration (hydrofoil, cylinder, bend, sphere, sphere+URF, across serial_regression.py, parallel_regression.py and hybrid_regression.py's copy of the hydrofoil case) - all re-verified to pass at the standard 1e-5 tolerance. inc_flatplate_pb's registered values needed no change. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…byA's A_p SetMomCoeff's A_p already has Vol/dt added to it (from PrepareImplicitIteration_impl) before SIMPLEC's A_p-Sum_A_nb correction runs. With TRANSIENT_TERM_REMOVAL_FACTOR left at its 0.0 default, that correction evaluates to roughly Vol/dt regardless of the momentum operator, so the pressure correction becomes vanishingly weak at low CFL - silently defeating the point of choosing SIMPLEC over SIMPLE. CConfig::SetPostprocessing now defaults the factor to 1.0 when KIND_PB_ITER= SIMPLEC is selected and the option was not set explicitly (checked via OptionIsSet, the same pattern already used for OUTPUT_FILES's default), with a one-line warning naming the substitution so it is not silent. An explicit user value, including an explicit 0.0, is always respected. ComputeHbyA separately recomputed A_p = Jacobian(i,i)(1,1)/rho raw, rather than reading the value SetMomCoeff already produced - so once SIMPLEC or a nonzero removal factor is in play, the PISO correction and the pressure equation would use two different, inconsistent momentum coefficients. Replaced it with Vol/GetMomCoeff(iPoint), the same corrected A_p SetMomCoeff stores (which also picks up F1's strong-BC reconstruction for free). Hoisted the computation out of the iDim loop it was needlessly repeated in - it does not depend on iDim. Both changes are confirmed no-ops on the entire current pressure-based suite: none of the nine registered cases use KIND_PB_ITER= SIMPLEC, and ComputeHbyA's fix is bit-for-bit unchanged even at full Cauchy convergence on the cylinder case (checked to every printed digit at iteration 425), exactly as both cards predicted for a SIMPLE-only suite with the removal factor at its default. Exercising them together (SIMPLEC, removal factor left unset so the new default applies) on a scratch cylinder variant surfaced a new, undiagnosed finding instead of the expected improvement: KIND_PB_ITER= SIMPLEC itself diverges to unphysical values by inner iteration 2, at both the case's registered CFL and at CFL=1 - ruling out the CFL-sensitivity explanation that resolved every other stability issue found so far in this pass. With F3, F4 and F5 all in place, this is a separate, undiagnosed defect. No SIMPLEC test case was added to the registered suite as a result of this finding. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…ting it failed ApplyPressureVelocityCorrection's edge-loop mass-flux correction read poisson_nodes->GetGradient_Primitive(iPoint,0,iDim), which always returns 0.0 - CScalarVariable never overrides the 3-argument form, only CFlowVariable does - silently degenerating CorrectPressureGradient to the compact two-point gradient with no non-orthogonal correction, despite the surrounding code (and the OpenFOAM/uFVM reference it cites) clearly intending one. Traced the actual mechanism before assuming a one-line accessor fix would activate it correctly: poisson_nodes->GetGradient() (the non-primitive accessor) is not stuck at zero the way GetGradient_Primitive is - it is correctly populated by CPoissonSolver::Postprocessing after each PISO correction's linear solve, and is already used correctly a few lines above this bug for the interior per-point velocity correction. Only the edge-loop mass-flux correction read the wrong accessor. Tried activating it (GetGradient instead of GetGradient_Primitive), reasoning that since the interior correction already legitimately uses this same real gradient with no apparent issue, extending it to the edge loop should be safe. It was not: the 3D sphere case diverges from inner iteration 1 regardless of CFL, with a sign-flipping, magnitude-growing CD trajectory - the same oscillation signature as the unrelated broken-SIMPLEC finding. The cylinder case also measurably degraded. The likely mechanism is that a node-averaged gradient at a face is exactly the quantity Rhie-Chow interpolation exists to avoid using directly in a mass flux - averaging it back in reintroduces checkerboard-style pressure-velocity decoupling. Reverted, and landed the other option the underlying finding offers instead: made the always-zero behavior explicit (GradPressure_avg[iDim] = 0.0;) rather than reachable only by accident through an unrelated base-class stub, with a comment recording both why it must stay zero and that a future contributor who notices GetGradient_Primitive "looks like a bug" would silently reintroduce the exact instability just found here. Confirmed bit-for-bit unchanged behavior on all four registered pressure-based cases - this changes documentation and defensive intent, not behavior. Left CAvgGrad_Heat's correct=true construction in CDriver.cpp untouched: it is shared with CHeatSolver, and though the same underlying Poisson equation mechanism makes its own deferred correction structurally always-zero too (Preprocessing resets p'=0 before recomputing its gradient, every PISO correction), that only costs wasted cycles, and touching shared numerics construction for a performance-only concern was judged out of scope here. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The X-0 guard against KIND_INCOMP_SYSTEM= PRESSURE_BASED + TIME_DOMAIN= YES cited "the Rhie-Chow mass flux interpolation has no transient term" - a restatement of the original review finding, not something independently verified against this guard's own presence. Removed the guard temporarily and actually ran incomp_pb_cylinder.cfg with TIME_DOMAIN= YES and DUAL_TIME_STEPPING-2ND_ORDER to see what happens: it does not crash, it runs to completion, but produces physically meaningless output (CD swinging between roughly -9.8 million and +7.6 thousand across 5 time steps while rms[P] sits steady around -16.6) - a silently wrong answer, which is worse than a crash and makes the guard more important, not less. Traced the cause: CPBFluidIteration::Iterate never calls SetDualTime_Solver or SetDualTime_Geometry for any solver, and never inspects TimeMarching or Time_Domain at all, unlike the standard CFluidIteration::Iterate which drives that bookkeeping for the density-based path. So the gap isn't one missing term in Rhie-Chow, it's that the entire PB outer-iteration driver has no time-history awareness whatsoever - each "time step" just partially re-converges the same steady problem from wherever the previous one left off, uncoupled from dt, matching the wild CD swings observed. Updated the guard message to state the actual, verified reason instead of the original restated hypothesis. Not attempting to implement dual-time support itself - that's a real feature addition (wiring the standard dual-time machinery into a driver that was never built with it, and deciding what "previous time step" means for a correction variable like p' that resets every PISO correction) rather than a bounded fix. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The previous message claimed CPBFluidIteration has no time-history bookkeeping at all. That was wrong: CPBFluidIteration only overrides Iterate() - Update() (which shifts Solution_time_n/Solution_time_n1) is inherited unchanged from CFluidIteration and does run, and SetResidual_DualTime is invoked transparently by the shared CIntegration::Space_Integration used by every MultiGrid_Iteration/SingleGrid_Iteration call regardless of which iteration driver calls it. So the physical dual-time source term for the momentum equation is already being added automatically. With that corrected, re-investigated what actually breaks: removed the guard again and ran incomp_pb_cylinder.cfg with TIME_DOMAIN= YES. It does not diverge - rms[U]/rms[V] converge cleanly below -10 within each time step - but converges to a self-consistent, physically wrong state. Dumped the raw converged fields directly: the pressure field swings to roughly +-965 and velocity reaches 1.88, despite an initial condition and far-field value of 0.000008 (six orders of magnitude smaller) with no forcing anywhere - a fluid at rest with all boundaries also at rest should stay at rest. Ruled out a simple physical-timestep-scaling explanation by rerunning at TIME_STEP=100 (10,000x larger): result did not converge toward the correct answer and even flipped sign between time steps. Root cause not established - likely candidates are an interaction between the dual-time Jacobian contribution and SetMomCoeff's A_p (which Coeff_Mom and the Poisson diffusion coefficient both derive from), or the literature-documented need for an explicit unsteady Rhie-Chow correction term the original review finding described - but neither confirmed. Updated the guard message to state only what was actually established rather than either the original restated hypothesis or this session's own earlier, incorrect "no bookkeeping" claim. Not attempting an implementation: this needs systematic term-by-term diagnosis and validation against a temporal-refinement study, not further black-box trial and error. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
- Gate the PB config-validation guard block on Kind_Regime == INCOMPRESSIBLE, not just the raw KIND_INCOMP_SYSTEM value. That option is read regardless of solver family, so a compressible or SU2_DEF config carrying a leftover PRESSURE_BASED line would otherwise hard-error on an unrelated MGLEVEL/TIME_DOMAIN/adjoint check. - Inline the substance of the TIME_DOMAIN guard's explanation instead of pointing at PB_SOLVER_PLAN.md, which is a local untracked planning file and not part of this branch's committed tree - anyone else checking out pr2812 would have hit a dangling reference. - Fix CIncNSSolver::Preprocessing's primitive-gradient switch to also handle plain NUM_METHOD_GRAD= LEAST_SQUARES, matching the fix adc77f6 already applied to the inviscid path (CIncEulerSolver::Preprocessing / ComputeEdgeMassFluxesRhieChow). The viscous solver calls CommonPreprocessing directly rather than the Euler override, so it had its own copy of the same gap: with plain LEAST_SQUARES the pressure gradient silently stayed zero. Confirmed no registered incompressible Navier-Stokes regression case uses plain LEAST_SQUARES, so this doesn't move any existing test_vals. - Add explanatory comments to CPBConvection_Base::ComputeJacobian's continuity and enthalpy rows, which still carry proj_vel-derived terms that contradict the "mass flux is frozen, so d(m_f)/du = 0" reasoning 9592182 used to justify dropping the equivalent momentum-row term. Left the values as-is rather than changing them blind: the continuity row is provably harmless (PrepareImplicitIteration_impl deletes row 0 for the pressure-based solver before the linear system is assembled), and the enthalpy row/column only matter when INC_ENERGY_EQUATION= YES, which no currently-registered pressure-based regression case exercises - so there's no test coverage to validate a re-derivation against. - Disable (rather than delete) the tutorials.py lid_driven_cavity entry: su2code/Tutorials#86, which adds the referenced Inc_Lid_Driven_Cavity/incomp_pb_liddrivencavity.cfg and its mesh, has not landed (checked the local Tutorials clone - no matching branch or content). Leaving the entry active would fail CI outright rather than fail a regression check. Kept the recorded test_vals in a comment so whoever lands that Tutorials PR can re-enable this with one uncomment, not a re-derivation. Not changed, on reflection: - The RowDeleted second Jacobian-diagonal sweep in SetMomCoeff (flagged as an efficiency finding) does add real cost, but a correctness-safe fix would need CIncEulerVariable::strongBC extended to cover wall and outlet row-deletions too (currently only set for the far-field-inflow case), which is a real behavioral change I'm not making blind. Left as-is. - The suspected stale-halo-density issue in ComputeHbyA across PISO corrections doesn't hold up: CIncEulerVariable::SetPressure/SetVelocity only touch the Pressure/Velocity primitive entries, never Density, so density stays frozen at its "top of outer iteration" value for every point (local or halo) throughout the correction loop - there's no actual divergence to fix. Verified Common/src/CConfig.cpp, SU2_CFD/src/solvers/CIncNSSolver.cpp and SU2_CFD/src/numerics/flow/convection/pressure_based.cpp compile clean, the full SU2_CFD binary links, and incomp_pb_NACA0012.cfg (GREEN_GAUSS, unaffected by the LEAST_SQUARES fix) still runs and converges to Exit Success. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CPBConvection_Base::ComputeJacobian assigned nonzero d(EnergyFlux)/d(velocity) and d(MomentumFlux)/d(enthalpy) entries, even though Flux[nDim+1] = MassFlux * AdvectedEnthalpy (frozen MassFlux, AdvectedEnthalpy an average of h_i/h_j only) has zero dependence on velocity, and Flux[1+iDim] = MassFlux * AdvectedVelocity has zero dependence on enthalpy - by the exact same "frozen mass flux" argument the code already uses to justify the momentum block's own Jacobian a few lines above. The comment already sitting on this code flagged the issue but nobody had run a case to confirm it mattered. Confirmed empirically: built a scratch coupled-energy PB case (viscous cylinder, INC_ENERGY_EQUATION=YES, isothermal wall) and compared against the same physics on the density-based path. Before this fix, the PB path failed to converge (rms[U] stalled at -0.38, 1000+ non-physical points, temperature swinging 0K-5000K) regardless of CFL (tried 50 and 1, lower was worse, ruling out a stability-tuning explanation) while the density-based path converged cleanly to a sane field (288.0-288.5K). After zeroing the two spurious cross-coupling entries, the PB path's temperature field settles to 288.1-288.6K, matching the density-based reference almost exactly, and the non-physical-point warnings are gone. A separate, milder issue remains: even fixed, the case converges only to an oscillating rms~-1 to -1.6 plateau rather than the density-based path's clean convergence - plausibly a missing energy-specific relaxation factor (F2 only covers pressure/velocity), not confirmed as a code bug and not addressed here. No effect on any currently-registered PB regression case: all four already run with INC_ENERGY_EQUATION=NO, and the existing !energy branch already zeroes the entire row/column nDim+1 after ComputeJacobian runs, independent of this change. Re-verified bit-for-bit identical on all four (cylinder, NACA0012, bend, sphere) via the TestCase.py harness. Separately confirmed (no code change needed): weakly-coupled heat mode (WEAKLY_COUPLED_HEAT=YES) is already correctly wired up for PB, since CPBFluidIteration::Iterate calls the same shared CommonAuxiliarySolvers that dispatches HEAT_SOL as the density-based path. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
… density) Adds pb_poly_cylinder.cfg (isothermal-wall heated cylinder, polynomial viscosity/conductivity, INC_ENERGY_EQUATION=YES, INC_DENSITY_MODEL=VARIABLE) and registers it in both parallel_regression.py and serial_regression.py at test_iter=20, with separately-verified serial/parallel test_vals (they differ slightly due to partitioning). This is the first regression coverage for the coupled-energy Jacobian fix landed earlier (bfd4181) - previously that fix was only validated via scratch configs that never got checked in, so nothing would catch a regression there. It's a short trajectory guard, not a converged-solution check: rms[h] converges genuinely (monotonically, no oscillation) but very slowly on this case - reaching a residual like -10 would need on the order of 40,000 iterations, so this test only pins down the first 20 iterations' behavior, the same pattern used elsewhere in this suite for cases that are correct but slow to converge. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Three review findings, none of which change results: Comments described the investigation that produced the code rather than the mechanism a reader needs: what the entries used to be, what was tried, which case diverged and by how much. Cut to the mechanism. Two of them also named the (uncommitted, repo-external) planning document from source, which must leave no trace in the repository. The TIME_DOMAIN guard printed a fifteen-line writeup at the user, including residual values and candidate root causes. A user hitting it needs to know the combination is unsupported and that it fails silently rather than loudly; the rest belongs outside the source. SetMomCoeff and ComputeEdgeMassFluxesRhieChow are consecutive with nothing serial between them, but each opened its own parallel region. Merged into one. CompleteComms ends in SU2_OMP_SAFE_GLOBAL_ACCESS, which barriers on both sides, so the write-then-read between the two calls is still ordered. Verified unchanged: all pressure-based cases in the parallel, serial and hybrid lists, the last at the -t 2 the hybrid suite actually runs. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The momentum equation is not assembled at a strong velocity BC: DeleteValsRowi zeroes the row and writes 1.0 on the diagonal, so there is no A_p to read at a wall, inlet or inflow far-field point. SetMomCoeff reconstructed one by averaging the raw coefficients of the point's neighbours, staging two nPoint arrays and scanning every point's Jacobian row against every neighbour to decide which points needed it. That extrapolates from an operator the boundary condition overwrote, and it is not what the coefficient is used for. Every consumer that needs a real value is edge-based - the Poisson diffusion coefficient and the two Rhie-Chow terms - so an edge touching such a point can simply use the coefficient of its other node. The two point-based consumers need no value at all: the velocity correction is overwritten in the boundary loop for each of those marker types, and HbyA sums A_nb over a deleted row, which is identically zero whatever it is divided by. So the reconstruction is dropped. Points under a strong velocity BC are identified by the strongBC flag, which already existed and was already allocated for this solver but was only set at inflow far-field; it is now also set where the wall and inlet conditions delete their rows, guarded to the pressure-based path because the flag is not allocated otherwise. They store a finite placeholder, which also avoids dividing by zero in the SIMPLEC correction at a transient removal factor of 1. RawMomCoeff, RowDeleted and the whole first pass are gone. Evidence, on the laminar bend (walls, inlet, outlet and symmetry all monitored), each run to convergence on the same mesh: density-based reference CD = 1.7069 (iter 1316) neighbour averaging CD = 0.3990 (iter 110) -77% neighbour of the edge CD = 1.9425 (iter 113) +14% Results move on every pressure-based case, most on the bend, by a few decimals elsewhere; reference values are updated in all three lists, the hybrid ones at the -t 2 that suite runs. Known gap, unchanged by this commit: strongBC is reset over nPointDomain only, so a halo point always reads false and an edge whose off-rank node is a wall will not see it. The averaging pass had the same blind spot, it skipped off-rank neighbours. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The OpenMP corrections to this solver - the Poisson boundary loop running on every thread, the reducer-strategy races, the stages that ran outside any parallel region - were covered by a single hybrid case, the inviscid hydrofoil. It exercises neither the viscous path, nor 3D, nor most of the marker types in the Poisson source loop. Adds the laminar cylinder (viscous, wall markers), the sphere (the only 3D case, and the one where the alpha_p reduction sums a different set of Jacobian diagonals per thread) and the heated cylinder (coupled energy). Values recorded with -t 2, the thread count the hybrid suite runs at; they differ from the MPI ones, so recording them at any other thread count would have failed CI. Also drops the poly cylinder and bend duplicates from the serial list. Running the same pressure-based case serially and in parallel buys a partitioning difference, which is worth having for the cases that were already in both, but was not a reason to double every case added since. The split is now: parallel runs the case, hybrid exercises the threaded path. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The edge loop is colored, which is what guarantees that no two edges handled concurrently share a point, so writing the source straight into LinSysRes at iPoint and jPoint is already safe. It is the same guarantee the viscous residual assembly relies on in this solver. Staging the per-edge source into an nEdge array and scattering it afterwards in a second loop over every point and its edges therefore bought nothing, at the cost of the array and of visiting every edge twice. Removed, along with EdgeSourceFlux. The mass flux and the HbyA correction are summed and added once rather than added separately, which is the only difference left from the original loop. The other half of the commit that introduced this, adding the work-sharing construct to the boundary marker loops that were running their full vertex range on every thread, is a real fix and is untouched. No change to any reference value in the three lists. 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
1 participant
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.
No description provided.