Lid driven cavity flow tutorial - #86

- 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>

## Proposed Changes 

The current work (part of GSoC) provides a working version of a
pressure-based algorithm for the incompressible flow solver as an
alternative to the existing Density-based solver. Below, the reader may
find the algorithm which has been implemented, as well as the current
progress of the code and challenges. All the way at the bottom one can
find performance comparisons between the DB and PB solvers for some test
cases.

### Algorithm

A lot of versions of pressure-based algorithms exist, and many different
versions can be implemented. Here, we opt for versions of the original
SIMPLE/PISO algorithm, it is briefly defined here for clarity.

First, the momentum equations are solved, starting from the previous
time step's velocity $\vec{u}^{(0)}$, pressure $p^{(0)}$, and face
velocity $\vec{u}_f^{(0)}$. The resulting momentum is the predicted
momentum, here its discretized form is shown, as its coefficients are
used in the subsequent equations

$A_p(\rho\vec{u})_p^{(1)}+\sum_n A_n (\rho\vec{u})_n^{(1)}=-V\nabla
p^{(0)}+S_m$

The subsequent momentum is not necessarily incompressible, the pressure
correction equation can be derived by rewriting it as follows, using a
term often called H by A

$(\rho\vec{u})_p^{(1)}=-\frac{\sum_n A_n
(\rho\vec{u})_n^{(1)}}{A_p}-\frac{V}{A_p}\nabla
p+S_m=\frac{H((\rho\vec{u})^{(1)})}{A}-\frac{V}{A_p}\nabla p^{(0)}+S_m$


$(\rho\vec{u})_p^{(2)}=\frac{H((\rho\vec{u})^{(1)})}{A}-\frac{V}{A_p}\nabla
p^{(1)}+S_m$

Note how a simplification is used here where the HbyA term is neglected.
Subtracting these two equations yields the first pressure correction
equation for $p'$ as

$\nabla \cdot \left( \frac{V}{A_p}\nabla p'\right)=\nabla \cdot \left(
\rho \vec{u}^{(1)}\right)$

Make note that for the divergence here, we require the face mass fluxes,
which are computed using Rhie-Chow interpolation to avoid odd-even
decoupling. After the equation is solved, using the pressure correction
$p'$, the pressure and momentum are corrected according to

$p^{(1)} = p^{(0)} + p' ,\quad (\rho u)^{(2)} = (\rho u)^{(1)}+(\rho
u)'.\quad (\rho u)'= -\frac{V}{A_p}\nabla p'$

So far, this is equal to a pseudo-transient version of the SIMPLE
algorithm. This algorithm however suffers from a very tight stability
condition on the time-step size. Therefore, multiple pressure
corrections can be applied, which for two corrections is originally
called the PISO algorithm.

The second pressure correction does not neglect the HbyA term, which
then results in the equation

$\nabla \cdot \left( \frac{V}{A_p}\nabla p'\right)=\nabla \cdot \left(
\rho \vec{u}^{(2)}\right)+\nabla\cdot\left(\frac{H((\rho
u)')}{A}\right)$

And the new correction equations are defined as

$p^{(2)} = p^{(1)} + p' ,\quad (\rho u)^{(3)} = (\rho u)^{(2)}+(\rho
u)'.\quad (\rho u)'= \frac{H((\rho u)')}{A}-\frac{V}{A_p}\nabla p'$

Note that HbyA here uses the previous velocity correction and is thus
the same quantity as the one used in the second pressure correction
equation. Later pressure correction equations follow analogously.


### Progress:

- Pressure-based solver added as alternative to density-based solver for
the incompressible flow equations.
- The pressure-based solver has only been tested for constant density
cases for basic Navier-Stokes and Euler flow.
- The pressure-based solver is implemented based on a pseudo
time-stepping approach to remain consistent with the other solvers in
SU2. The pressure-based solver is currently set to the SIMPLE algorithm
by default, with options for SIMPLEC and PISO available.
- The Poisson solver is a major bottleneck in the computation speed. To
account for this, an option was added to use a different linear solver
and preconditioner for the poisson solver.
- For details on the implementation of the algorithm and the
responsibility distribution please see the file CPBFluidIteration.cpp.
- A tutorial showcasing the different options of the PB solver for the
lid driven cavity flow problem has been added, please see
su2code/Tutorials#86 and su2code/su2code.github.io#218
- Regression tests have been  added.

### Future work:
#### Performance:

- The Poisson solver can sometimes struggle a lot due to high Reynolds
numbers and fine meshes, and thus require a ridiculous number of
iterations to converge reasonably. Possible fixes include adding
multigrid support or a DIC preconditioner (far less efficient).
Multigrid support is tricky as SU2 currently only considers multigrid
for the main (flow) solver and not for auxiliary solvers.

- Convergence issues with RANS (SA and SST) on fine meshes with high
Reynolds numbers. Tests have shown that cases such as flow over a flat
plate converges fine. However, external aerodynamic cases such as the
naca0012 RANS test case do not converge well at all. The convergence
does slightly improve when we switch out the mesh for a more uniform
unstructured mesh without large aspect ratio cells in the wake of the
airfoil, although this only slightly helps. The flat plate turbulence
test case also uses large aspect ratio cells so this is not the sole
issue. The RANS solver also often requires many iterations of the
Poisson solver to converge reasonably, this is however not the reason
for the lack of convergence.

- Periodic boundary conditions have not been implemented/tested at all
as of yet.
- Any code related to adjoints has not been considered at all either.

#### Code:

- Parallelization with OMP gives wrong results due to a unknown issue in
the Poisson solver, MPI however does work as expected.
- The restart solution currently does not write the edge mass fluxes
which are used by the pressure-based solver. Therefore restarted
solutions have to estimate these mass fluxes based on the average of the
nodal solutions. This results in the restarted solutions starting from
slightly different residuals. A fix for this issue can be to let the
edge mass fluxes be stored in the restart file.
- Support for the energy equation is added trough both the enthalpy
equation as well as trough the weakly coupled heat solver. This support
has however not yet been verified trough actual test cases.
- The numerics class for the convective residuals write the continuity
parts of the flux and jacobians. This will allow for a coupled solver in
the future, but as it stands now this code is _not_ used.

### TODO list 

- Fix issues mentioned above (left for future work)

## Related Work

This work is based on earlier attempts by Nitish Anand (2024) and Akshay
Koodly (2021), see feature branches feature_PBFlow_V8 and
feature_Pressure_based respectively. Also see PR #2210

## PR Checklist

- [x] I am submitting my contribution to the develop branch.
- [x] My contribution generates no new compiler warnings (try with
--warnlevel=3 when using meson).
- [x] My contribution is commented and consistent with SU2 style
(https://su2code.github.io/docs_v7/Style-Guide/).
- [x] I used the pre-commit hook to prevent dirty commits and used
`pre-commit run --all` to format old commits.
- [x] I have added a test case that demonstrates my contribution, if
necessary.
- [x] I have updated appropriate documentation (Tutorials, Docs Page,
config_template.cpp), if necessary.

## Result showcase 

### Inviscid Hydrofoil

<img
src="https://github.com/user-attachments/assets/2178b449-557d-4ec5-89b6-0d49ace6b612"
width="600">

*Convergence history of the inviscid flow around a hydrofoil at a 5
degree aoa.*

<img
src="https://github.com/user-attachments/assets/cd28de03-1c0a-4927-a5a3-567c454eb5c2"
width="600">

*The pressure coefficient along the surface of the hydrofoil at a 5
degree aoa and the corresponding lift coefficients, X-FOIL predicts
C_L=0.6.*

### Lid Driven Cavity 

<img
src="https://github.com/user-attachments/assets/0433d13c-b76d-4b9d-ac2c-615539877c9e"
width="600">

*Convergence history of the lid driven cavity problem, note that CFL=60
is the highest stable CFL for the PB solver, whereas the DB solver does
not have this CFL related stability issue.*

### Flatplate RANS

<img
src="https://github.com/user-attachments/assets/b8db128d-e919-46a5-9b9d-aa17c1b94219"
width="600">

*The skin friction coefficient for turbulent flow over a (rough) flat
plate with SA.*

---------

Co-authored-by: Pedro Gomes <pcarruscag@gmail.com>
Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>