Cache Least Squares matrix for gradient computation

[tbellosta]

Proposed Changes

This PR adds an opt-in option, CACHE_LSQ_METRICS (default NO), caches the factorized metric terms of the (weighted or unweighted) least-squares gradients. The LSQ normal matrix A depends only on the node coordinates and on the weighting, yet it is currently re-assembled and re-factorized on every gradient evaluation, every nonlinear iteration. With the option enabled:

• On the first evaluation the normal matrix is assembled and factorized as usual, and its inverse S = A^-1 (6 unique entries per node in 3D, 3 in 2D) is stored in CGeometry.
• Every subsequent evaluation only accumulates the right-hand side b = sum_k wd(u_k - u_i) in a loop over the edges (using the edge coloring; the contribution of an edge is identical for both of its end points, so each edge is visited once) followed by one small symmetric matrix-vector product S*b per node.
• The cache is stored per grid level with one slot per weighting, shared by all solvers (flow primitive and reconstruction gradients, turbulence/species/scalar solution gradients, auxiliary variables). The common practice of unweighted LSQ for MUSCL reconstruction plus inverse-distance-weighted LSQ for viscous gradients is fully supported: both slots are built once and reused; e.g. the turbulence solver directly reuses the metrics built by the flow solver.
• Moving/deforming grids are supported: the cache is invalidated in SetControlVolume(..., UPDATE) (the common point of all mesh motion/deformation paths, on fine and coarse MG levels) and rebuilt on the next evaluation, i.e. once per mesh update, so the savings scale with the number of inner iterations per time step.
• The option is automatically disabled (with a warning) for periodic boundaries (the periodic LSQ communication fuses the matrix and RHS accumulations; an RHS-only exchange would be needed) and for the discrete adjoint (the coordinate dependence of the metrics must remain on the tape, just a guess I am not an expert on the discrete adjoint implementation).

Saved operations (per node per gradient evaluation, 3D, k edge-neighbors, N variables)

Eliminated from every evaluation after the first:

• normal-matrix assembly: ~21k flops (~300 flops per node on tetrahedra, k=14),
• per-node Cholesky factorization + inversion: ~35 flops and 3 sqrt (no sqrt/div remain in the iteration path),
• half of the geometric edge work (distance vectors and solution differences are computed once per edge instead of once per node-neighbor pair).

What remains is the RHS accumulation (~4N flops per node per edge) and one S*b product (18N flops per node). For the N=6 compressible primitives this roughly halves the gradient-kernel flops; for scalar solvers (N=1-2, e.g. SA) where the assembly dominates it is a ~4-5x reduction. Measured end-to-end (serial, linear-solver work held fixed): ~3.3% of total wall time on an inviscid ONERA M6 (57.5k nodes, tets, one gradient set per iteration) and ~6% on a 2D RANS-SA case (three LSQ gradient sets per iteration). Cases dominated by the linear solver will see proportionally less.

Memory footprint

The cache adds nDim*(nDim+1)/2 su2doubles per node per weighting used (48 B/node per weighting in 3D, 24 B/node in 2D; at most two weightings), allocated lazily per grid level and shared by all solvers. For reference, every solver already allocates an nDim*nDim Rmatrix scratch (72 B/node in 3D), so for a typical RANS case the addition is small compared to the existing gradient machinery and negligible compared to overall solver memory.

Validation

With the option enabled, results are identical to develop (to output precision) for:

• inviscid NACA0012, unweighted LSQ reconstruction, 200 iterations (serial and 2 MPI ranks),
• RANS-SA NACA0012 with NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES and NUM_METHOD_GRAD_RECON= LEAST_SQUARES (all residuals including the turbulence variable),
• an unsteady pitching NACA0012 with GRID_MOVEMENT= RIGID_MOTION (dual time stepping), exercising the invalidation/rebuild path,
• hybrid-parallel runs (OpenMP edge-colored loop, bit-reproducible run-to-run; falls back to a race-free node loop if the natural coloring is in use).

config_template.cfg documents the new option. Default behavior (CACHE_LSQ_METRICS= NO) is bit-identical to current develop.

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

[pcarruscag]

Nice. Make it the default, compute in geometry preprocessing when necessary, no opt-in

[pcarruscag]

The failed tests seem to be just due to floating-point ordering. Is the PR ready for review?

[tbellosta]

The failed tests seem to be just due to floating-point ordering. Is the PR ready for review?

Kind of. I wanted to run some more detailed measurements on larger runs and add the caching for periodic boundaries as well. I will have some spare time next week to do that. I also have to run the pre-commit hook

[tbellosta]

I was finally able to implement the caching for period boundaries and run some benchmark cases.

Periodic boundaries

The periodic contributions to the least-squares normal matrix are exchanged at the geometry level (a metric-only periodic communication added to the geometry, independent of any solver), so the cache is built once in CDriver::InitializeGeometry for periodic and non-periodic grids alike. The solver-level periodic communication of the LSQ matrix is skipped when caching is active, only the RHS still goes through the usual solver periodic comm. Results on the periodic regression cases (poiseuille profiles, channel, periodic test cases of the parallel/hybrid suites) are identical to develop.

Edge-coloring fix

An earlier version of this branch built the edge coloring eagerly in the cache setup,
before the solver constructors, which changed the solvers' reducer-strategy decision on
multi-rank runs and gave slightly different results versus develop with np > 1. The
coloring is now left to the solvers as before, the cached gradient path only builds it
itself when no solver did. With this fix the new binary matches develop on
1 to 4 ranks, 1 thread per rank, and in hybrid mode.

Benchmark

Full regression suites (parallel, hybrid, serial) run with both the develop base
(b12aa0c) and this branch on the same node (AMD EPYC 9654, gcc 11.2, -march=native, OpenMPI 4.1.6). I had to exclude NEMO cases since I wan not able to run them even with the upstream develop.

	develop	this PR	change
Sum of compute time, 332 cases	805.2 s	784.6 s	-2.6 %
Sum of wall time	1345.1 s	1288.1 s	-4.2 %
Sum of preprocessing time	129.7 s	126.8 s	-2.3 %
Compute, parallel suite (151 cases)			-2.3 %
Compute, hybrid suite (72 cases)			-2.3 %
Compute, serial suite (109 cases)			-2.9 %

Of the 147 cases with more than 1 s of compute time, 124 are faster, 19 are within the
±1 % noise band and 4 are slower (DG-FEM cases that do not use the LSQ path, so noise).
Largest gains: cfd_flamelet_ch4_unsteady -11 %, dyn_fsi -10 %, fsi2d -8 %,
fem_ns_sphere -8 %, cfd_flamelet_ch4 -7 %.

Note on regression references

The cached path accumulates the RHS in an edge loop instead of in a loop through point neighbours, which changes the roundoff of the gradients. This will make some regression cases fail. Those references will need to be updated once CI runs on this branch, unless we prefer to keep neighbour-order accumulation in the cached path at some cost in speed.

IMO this is now ready for review.