Multigrid projection integration

CHANGELOG.md

@@ -9,6 +9,25 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- **Multigrid wired into the projection method** — new
+  `ns_solver_params_t.pressure_solver` field (`ns_pressure_solver_t`; 0 =
+  existing CG behavior) selects the pressure Poisson solve of the scalar
+  `projection` solver: `NS_PRESSURE_SOLVER_MULTIGRID` (multigrid V-cycles,
+  RHS interior mean subtracted for Neumann compatibility) or
+  `NS_PRESSURE_SOLVER_PCG_MG` (MG-preconditioned CG). Non-2^k+1 grids and
+  non-scalar projection backends reject the selection with
+  `CFD_ERROR_UNSUPPORTED` at init
+  (`lib/src/solvers/navier_stokes/cpu/solver_projection.c`,
+  `lib/src/api/solver_registry.c`,
+  `tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c`).
+- **Multigrid-preconditioned CG** — `POISSON_PRECOND_MULTIGRID` runs one
+  symmetric V(2,2) weighted-Jacobi multigrid cycle in Dirichlet mode per
+  preconditioner apply (scalar CG only; other CG/GMRES backends reject it
+  with `CFD_ERROR_UNSUPPORTED`). Grid-size-independent convergence: 5 CG
+  iterations at 33²–129² (tol 1e-8) vs 50–170 unpreconditioned. Exposed as
+  the `POISSON_SOLVER_PCG_MG_SCALAR` convenience preset
+  (`lib/src/solvers/linear/cpu/linear_solver_cg.c`,
+  `tests/math/test_mg_pcg_convergence.c`).
 - **Geometric multigrid Poisson solver** (scalar backend) — V/W/F(FMG) cycles with
   Red-Black Gauss-Seidel or weighted-Jacobi smoothers, full-weighting restriction and
   bilinear/trilinear prolongation, 2D/3D. Two BC modes: Neumann zero-gradient (default,

CMakeLists.txt

@@ -253,6 +253,7 @@ if(BUILD_TESTS)
     add_executable(test_solver_rk2 tests/solvers/navier_stokes/cpu/test_solver_rk2.c)
     add_executable(test_solver_rk4 tests/solvers/navier_stokes/cpu/test_solver_rk4.c)
     add_executable(test_solver_projection tests/solvers/navier_stokes/cpu/test_solver_projection.c)
+    add_executable(test_projection_pressure_solver tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c)
     add_executable(test_ns_solver_3d tests/solvers/navier_stokes/cpu/test_ns_solver_3d.c)
 
     # Navier-Stokes solver tests (AVX2)
@@ -348,6 +349,7 @@ if(BUILD_TESTS)
     add_executable(test_divergence_free tests/math/test_divergence_free.c)
     add_executable(test_mms tests/math/test_mms.c)
     add_executable(test_pcg_convergence tests/math/test_pcg_convergence.c)
+    add_executable(test_mg_pcg_convergence tests/math/test_mg_pcg_convergence.c)
     add_executable(test_residual_computation tests/math/test_residual_computation.c)
     add_executable(test_cg_scaling tests/math/test_cg_scaling.c)
     add_executable(test_nonuniform_grid tests/math/test_nonuniform_grid.c)
@@ -407,6 +409,7 @@ if(BUILD_TESTS)
     target_link_libraries(test_solver_rk2 PRIVATE CFD::Library unity)
     target_link_libraries(test_solver_rk4 PRIVATE CFD::Library unity)
     target_link_libraries(test_solver_projection PRIVATE CFD::Library unity)
+    target_link_libraries(test_projection_pressure_solver PRIVATE CFD::Library unity)
     target_link_libraries(test_ns_solver_3d PRIVATE CFD::Library unity)
     target_link_libraries(test_solver_explicit_euler_avx2 PRIVATE CFD::Library unity)
     target_link_libraries(test_solver_projection_avx2 PRIVATE CFD::Library unity)
@@ -478,6 +481,7 @@ if(BUILD_TESTS)
     target_link_libraries(test_divergence_free PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_mms PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_pcg_convergence PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
+    target_link_libraries(test_mg_pcg_convergence PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_residual_computation PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_cg_scaling PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_nonuniform_grid PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
@@ -521,6 +525,7 @@ if(BUILD_TESTS)
     add_test(NAME SolverRK2Test COMMAND test_solver_rk2)
     add_test(NAME SolverRK4Test COMMAND test_solver_rk4)
     add_test(NAME SolverProjectionTest COMMAND test_solver_projection)
+    add_test(NAME ProjectionPressureSolverTest COMMAND test_projection_pressure_solver)
     add_test(NAME NSSolver3DTest COMMAND test_ns_solver_3d)
     add_test(NAME SolverExplicitEulerAvx2Test COMMAND test_solver_explicit_euler_avx2)
     add_test(NAME SolverProjectionAvx2Test COMMAND test_solver_projection_avx2)
@@ -642,6 +647,7 @@ if(BUILD_TESTS)
     add_test(NAME DivergenceFreeTest COMMAND test_divergence_free)
     add_test(NAME MmsTest COMMAND test_mms)
     add_test(NAME PCGConvergenceTest COMMAND test_pcg_convergence)
+    add_test(NAME MgPcgConvergenceTest COMMAND test_mg_pcg_convergence)
     add_test(NAME ResidualComputationTest COMMAND test_residual_computation)
     add_test(NAME CgScalingTest COMMAND test_cg_scaling)
     add_test(NAME NonuniformGridTest COMMAND test_nonuniform_grid)

README.md

@@ -12,7 +12,7 @@ A production-grade computational fluid dynamics (CFD) library in C for solving 2
 
 - 🚀 **Multiple Backends**: CPU (scalar), SIMD (AVX2/NEON), OpenMP, CUDA
 - 🔧 **Pluggable Solvers**: Explicit Euler, RK2 (Heun), Projection Method (Chorin's algorithm)
-- 📊 **Linear Solvers**: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB
+- 📊 **Linear Solvers**: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB, GMRES(m), geometric multigrid (V/W/F cycles, also as CG preconditioner and projection pressure solver via `params.pressure_solver`)
 - 🌡️ **Heat Transfer**: Energy equation (advection–diffusion) + Boussinesq buoyancy + thermal BCs
 - 🌀 **Turbulence (RANS)**: Standard k-ε and Spalart-Allmaras with log-law wall functions; validated vs turbulent channel flow at Re_τ = 395
 - 🎯 **Validated**: Ghia lid-driven cavity, Taylor-Green vortex, Poiseuille flow, natural convection, turbulent channel flow

ROADMAP.md

@@ -76,7 +76,9 @@ Genuine constraints to be aware of (not backlog items):
   Windows/macOS handle automatically. Future option: weak symbols or a plugin architecture.
 - **Multigrid grid-dimension constraint** — the geometric multigrid solver requires
   2^k+1 points per active dimension (e.g. 33, 65, 129) for exact coarsening; other sizes
-  return `CFD_ERROR_INVALID`. Use CG for arbitrary grid sizes.
+  return `CFD_ERROR_INVALID`. Use CG for arbitrary grid sizes. On conforming grids,
+  MG-preconditioned CG (`POISSON_PRECOND_MULTIGRID`) gives grid-size-independent
+  iteration counts on uniform grids where Jacobi PCG gives no benefit.
 
 > The SOR optimal-ω and Jacobi spectral-radius formulas (ρ = cos(πh),
 > ω = 2/(1 + sin(πh))) assume Dirichlet BCs; with Neumann BCs optimal ω is typically lower
@@ -130,8 +132,10 @@ non-symmetric operators arrive, e.g. implicit advection-diffusion in §1.5).
 - [x] Geometric multigrid — scalar backend; V/W/F(FMG) cycles, Red-Black GS or weighted
       Jacobi smoothers, full-weighting restriction (Neumann-folded at boundaries) +
       bilinear/trilinear prolongation, Neumann (default) and Dirichlet BC modes, 2D/3D.
-      Grid dims must be 2^k+1. SIMD/OMP/GPU variants and MG-preconditioned CG deferred
-      (see `.claude/specs/multigrid-projection-integration.md` for projection wiring).
+      Grid dims must be 2^k+1. Wired into the projection method
+      (`ns_solver_params_t.pressure_solver`, scalar backend) and available as a CG
+      preconditioner (`POISSON_PRECOND_MULTIGRID`, scalar CG). SIMD/OMP/GPU variants
+      deferred (see `.claude/specs/multigrid-projection-integration.md`).
 - [ ] Algebraic multigrid (AMG) — solver and preconditioner (for CG/GMRES/BiCGSTAB)
 - [x] GPU plain SOR (Block SOR: per-thread tile sweep, red-black tile coloring, in-place; closes the matrix)
 
@@ -209,7 +213,10 @@ SIMD Poisson integration is done; current ~1.3–1.5× speedup is Amdahl-limited
 (parallelizable fraction ~80%). Remaining optimization work:
 
 - [ ] Increase `POISSON_MAX_ITER` or implement adaptive tolerance
-- [ ] Optional multigrid preconditioner for faster convergence (see §1.2 multigrid)
+- [x] Optional multigrid preconditioner for faster convergence (see §1.2 multigrid) —
+      done for scalar CG (`POISSON_PRECOND_MULTIGRID`, symmetric V(2,2) Jacobi cycle) and
+      selectable in the scalar projection via `pressure_solver = NS_PRESSURE_SOLVER_PCG_MG`;
+      SIMD/OMP variants follow the multigrid backend work in §1.2
 - [ ] Red-Black omega parameter tuning
 - [ ] Profile to identify remaining bottlenecks
 - [ ] OpenMP+SIMD hybrid projection (OMP across rows, SIMD within rows); benchmark vs pure OMP

docs/guides/examples.md

@@ -681,12 +681,12 @@ CFD Platform Diagnostics
 - Error handling for unavailable solvers
 
 **Sections:**
-1. Method comparison (Jacobi, SOR, Red-Black SOR, CG, CG+Jacobi PC, BiCGSTAB) on scalar backend
+1. Method comparison (Jacobi, SOR, Red-Black SOR, CG, CG+Jacobi PC, CG+Multigrid PC, BiCGSTAB) on scalar backend
 2. Backend comparison (CG on Scalar, SIMD, OMP)
 3. Convenience API demo (`poisson_solve()`)
 4. Error handling (requesting multigrid on an unavailable backend — GPU)
 
-**Problem:** Solves ∇²p = -2π²sin(πx)sin(πy) on a 64×64 grid using the library's default homogeneous Neumann boundary conditions. The reported L2 error compares methods/backends against a common reference field — not against the Dirichlet analytical solution sin(πx)sin(πy), since BCs differ.
+**Problem:** Solves ∇²p = -2π²sin(πx)sin(πy) on a 65×65 grid (2^k+1, so the multigrid preconditioner can build its hierarchy) using the library's default homogeneous Neumann boundary conditions. The reported L2 error compares methods/backends against a common reference field — not against the Dirichlet analytical solution sin(πx)sin(πy), since BCs differ. Standalone multigrid is omitted: this RHS has a nonzero interior mean, which the true Neumann system cannot converge on (the same reason the stationary methods report max_iter).
 
 **Run:**
 ```bash
@@ -697,12 +697,18 @@ CFD Platform Diagnostics
 ```
 --- Method Comparison (Scalar Backend) ---
   Method                Iters     Residual    L2 Error      Time  Status
-  Jacobi                10001   res=8.3e+00  L2=4.8e+00  5548 ms  max_iter
-  CG                        1   res=5.1e-11  L2=1.1e-04     2 ms  converged
-  CG + Jacobi PC            1   res=5.9e-11  L2=1.1e-04     2 ms  converged
-  BiCGSTAB                  1   res=5.1e-11  L2=1.1e-04     2 ms  converged
+  Jacobi                10000   res=8.3e+00  L2=-1.0e+00   394 ms  max_iter
+  CG                        1   res=6.3e-11  L2=1.0e-04     0 ms  converged
+  CG + Jacobi PC            1   res=6.3e-11  L2=1.0e-04     0 ms  converged
+  CG + Multigrid PC         7   res=1.4e-06  L2=1.0e-04     2 ms  converged
+  BiCGSTAB                  1   res=6.3e-11  L2=1.0e-04     0 ms  converged
 ```
 
+> **Note:** `L2=-1.0e+00` is a "not computed" sentinel, not a real error value. The
+> example only computes the L2 error when the solve returns `CFD_SUCCESS`; on a
+> `max_iter` (non-converged) exit it prints `-1.0` instead. Jacobi does not converge
+> on this RHS within 10000 iterations, so its L2 error is reported as `-1.0`.
+
 ---
 
 ### 14. poiseuille_stretched_grid.c

docs/reference/api-reference.md

@@ -263,11 +263,26 @@ typedef struct {
     // Turbulence
     turbulence_model_t        turb_model;  // TURB_MODEL_NONE (default), K_EPSILON, SPALART_ALLMARAS
     ns_turbulence_bc_config_t turb_bc;    // Per-face turbulence BC types
+    // Pressure Poisson solver (projection solvers)
+    ns_pressure_solver_t pressure_solver;  // NS_PRESSURE_SOLVER_DEFAULT (0) = backend's CG
 } ns_solver_params_t;
 
 ns_solver_params_t ns_solver_params_default(void);
 ```
 
+`ns_pressure_solver_t` selects the pressure Poisson solve of the projection
+method (scalar `projection` solver only; the multigrid modes require 2^k+1
+grid points per active dimension, and other projection backends reject
+non-default values with `CFD_ERROR_UNSUPPORTED` at init):
+
+```c
+typedef enum {
+    NS_PRESSURE_SOLVER_DEFAULT = 0,   // Backend's default CG pressure solve
+    NS_PRESSURE_SOLVER_MULTIGRID = 1, // Geometric multigrid V-cycle solver
+    NS_PRESSURE_SOLVER_PCG_MG = 2,    // CG preconditioned by one MG V-cycle
+} ns_pressure_solver_t;
+```
+
 `turbulence_model_t` is defined in `cfd/solvers/navier_stokes_solver.h`:
 
 ```c
@@ -492,6 +507,13 @@ poisson_solver_params_t poisson_solver_params_default(void);
 > `poisson_solver_init` returns `CFD_ERROR_INVALID` otherwise. In the default
 > `MG_BC_NEUMANN` mode the solution is defined up to an additive constant.
 
+`poisson_precond_type_t` values: `POISSON_PRECOND_NONE` (0, default),
+`POISSON_PRECOND_JACOBI` (1), `POISSON_PRECOND_MULTIGRID` (2 — one MG V-cycle
+per apply; scalar CG only, 2^k+1 dims; other backends return
+`CFD_ERROR_UNSUPPORTED`). The convenience API exposes the MG-preconditioned CG
+as the `POISSON_SOLVER_PCG_MG_SCALAR` preset for
+`poisson_solve()`/`poisson_solve_3d()`.
+
 ### Poisson Statistics
 
 ```c
@@ -734,7 +756,7 @@ cfd_aligned_free(data);  // Correct
 #define NS_SOLVER_TYPE_PROJECTION_OPTIMIZED     "projection_optimized"
 #define NS_SOLVER_TYPE_PROJECTION_OMP           "projection_omp"
 #define NS_SOLVER_TYPE_RK2                      "rk2"
-#define NS_SOLVER_TYPE_PROJECTION_JACOBI_GPU    "projection_jacobi_gpu"
+#define NS_SOLVER_TYPE_PROJECTION_GPU           "projection_gpu"
 ```
 
 ## Version Information

docs/reference/solvers.md

@@ -76,7 +76,31 @@ Chorin's projection method - properly enforces incompressibility constraint.
 | `projection` | Scalar | Basic implementation |
 | `projection_optimized` | SIMD | SIMD-optimized (runtime detection: AVX2/NEON) |
 | `projection_omp` | OpenMP | Multi-threaded |
-| `projection_jacobi_gpu` | GPU | CUDA-accelerated (Jacobi iteration) |
+| `projection_gpu` | GPU | CUDA-accelerated (CG pressure solve) |
+
+**Pressure solver selection** (`ns_solver_params_t.pressure_solver`):
+
+Each backend pairs with a CG Poisson preset by default. On the scalar
+`projection` solver the pressure solve can be switched to geometric multigrid:
+
+| `pressure_solver` value | Pressure Poisson solve |
+|-------------------------|------------------------|
+| `NS_PRESSURE_SOLVER_DEFAULT` (0) | Backend's CG preset (existing behavior) |
+| `NS_PRESSURE_SOLVER_MULTIGRID` | Multigrid V-cycles (`POISSON_SOLVER_MG_SCALAR`) |
+| `NS_PRESSURE_SOLVER_PCG_MG` | CG + MG V-cycle preconditioner (`POISSON_SOLVER_PCG_MG_SCALAR`) |
+
+The MG modes require 2^k+1 grid points per active dimension (33, 65, 129, ...);
+`solver_init` returns `CFD_ERROR_UNSUPPORTED` otherwise. `projection_optimized`,
+`projection_omp`, and `projection_gpu` reject any non-default value with
+`CFD_ERROR_UNSUPPORTED` at init — multigrid has no SIMD/OMP/GPU backend yet and
+the library never falls back across backends silently.
+
+```c
+ns_solver_params_t params = ns_solver_params_default();
+params.pressure_solver = NS_PRESSURE_SOLVER_MULTIGRID;  // 2^k+1 grids only
+ns_solver_t* slv = cfd_solver_create(registry, NS_SOLVER_TYPE_PROJECTION);
+cfd_status_t status = solver_init(slv, grid, &params);   // UNSUPPORTED on 128x128
+```
 
 ## Linear Solvers (Poisson Equation)
 
@@ -196,6 +220,16 @@ poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_CG,
   - Simple, cheap per iteration
   - Variable coefficients: can reduce iterations
   - Constant coefficients (uniform grid): no benefit
+- **Multigrid** (`POISSON_PRECOND_MULTIGRID`): one geometric-multigrid V-cycle
+  per apply
+  - Grid-size-independent outer iteration count (measured 5 CG iterations at
+    33²–129², tol 1e-8, vs 50–170 for plain CG)
+  - Scalar CG backend only; other backends reject it with
+    `CFD_ERROR_UNSUPPORTED`
+  - Requires 2^k+1 grid points per active dimension (inherited from the
+    multigrid hierarchy); init fails with `CFD_ERROR_INVALID` otherwise
+  - The inner cycle is a symmetric V(2,2) with weighted-Jacobi smoothing in
+    Dirichlet mode, so the preconditioner is SPD as CG requires
 - **SSOR**: M = (D + L)D⁻¹(D + U) (future)
 - **ILU**: Incomplete LU factorization (future)
 
@@ -209,6 +243,19 @@ poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_CG,
 poisson_solver_init(solver, nx, ny, dx, dy, &params);  // Pass params with preconditioner
 ```
 
+**Multigrid-preconditioned CG (scalar backend, 2^k+1 dims):**
+
+```c
+poisson_solver_params_t params = poisson_solver_params_default();
+params.preconditioner = POISSON_PRECOND_MULTIGRID;
+
+poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_CG,
+                                                 POISSON_BACKEND_SCALAR);
+poisson_solver_init(solver, 65, 65, 1, dx, dy, 0.0, &params);
+```
+The convenience API exposes the same configuration as the
+`POISSON_SOLVER_PCG_MG_SCALAR` preset for `poisson_solve()`/`poisson_solve_3d()`.
+
 #### 6. BiCGSTAB
 
 **Algorithm:** Bi-Conjugate Gradient Stabilized for non-symmetric systems.
@@ -282,8 +329,10 @@ Smoothers: Red-Black Gauss-Seidel (default) or weighted Jacobi (ω=2/3).
 - Parameters: `mg_cycle`, `mg_smoother`, `mg_bc`, `mg_pre_smooth`/`mg_post_smooth`
   (default 2/2), `mg_coarse_max_iter` (default 50), `mg_max_levels` (0 = auto)
 
-**Backends:** scalar only. SIMD/OMP/GPU variants and use as a CG preconditioner
-are planned follow-ups.
+**Backends:** scalar only. Also available as a CG preconditioner
+(`POISSON_PRECOND_MULTIGRID`, scalar CG only — see §5) and as the projection
+method's pressure solver (`ns_solver_params_t.pressure_solver`). SIMD/OMP/GPU
+variants are planned follow-ups.
 
 **Usage:**
 ```c
@@ -307,6 +356,7 @@ poisson_solver_init(solver, 65, 65, 1, dx, dy, 0.0, &params);  // dims 2^k+1
 | Red-Black SOR | ~2000 | 8 | Parallelizable |
 | CG | ~80 | 5 | Best for large grids |
 | PCG (Jacobi) | ~80 | 5.5 | No benefit on uniform grid |
+| PCG (Multigrid) | ~5 | — | Grid-size-independent; needs 2^k+1 dims |
 | BiCGSTAB | ~40 | 4 | Fastest convergence |
 | Multigrid V(2,2) | ~8 cycles | — | O(N), grid-size-independent; needs 2^k+1 dims |
 
@@ -413,7 +463,7 @@ GPU-accelerated using CUDA:
 // Check if GPU should be used
 gpu_config_t config = gpu_config_default();
 if (gpu_should_use(&config, nx, ny, num_steps)) {
-    solver = cfd_solver_create(registry, "projection_jacobi_gpu");
+    solver = cfd_solver_create(registry, "projection_gpu");
 } else {
     solver = cfd_solver_create(registry, "projection_optimized");
 }
@@ -431,7 +481,7 @@ if (gpu_should_use(&config, nx, ny, num_steps)) {
 | projection | 19.0 | 1.0x | High |
 | projection_optimized | 5.3 | 3.6x | High |
 | projection_omp (8 cores) | 4.2 | 4.5x | High |
-| projection_jacobi_gpu | 8.4 | 0.45x† | High |
+| projection_gpu | 8.4 | 0.45x† | High |
 
 † GPU slower on small grids due to data transfer overhead
 
@@ -449,7 +499,7 @@ if (gpu_should_use(&config, nx, ny, num_steps)) {
 | Solver | Time (s) | Speedup |
 |--------|----------|---------|
 | projection_optimized | 824 | 1.0x |
-| projection_jacobi_gpu | 68 | 12.1x |
+| projection_gpu | 68 | 12.1x |
 
 ## Choosing a Solver
 
@@ -465,7 +515,7 @@ Need strict incompressibility enforcement?
 └─ Yes → Use Projection Method family
          ├─ Small grid (<100×100) → projection
          ├─ Medium grid (100-500) → projection_optimized or projection_omp
-         └─ Large grid (>500)     → projection_jacobi_gpu
+         └─ Large grid (>500)     → projection_gpu
 
 GPU available and grid >200×200?
 └─ Use CUDA variant for 10-50x speedup
@@ -480,7 +530,7 @@ GPU available and grid >200×200?
 
 **Production Simulations:**
 - `projection_optimized` or `projection_omp` (medium grids)
-- `projection_jacobi_gpu` (large grids)
+- `projection_gpu` (large grids)
 - Best accuracy and performance
 
 **Benchmarking:**