Add geometric multigrid Poisson solver (scalar backend)

CHANGELOG.md

@@ -9,12 +9,29 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- **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,
+  matching the other Poisson solvers; nullspace handled via Neumann-folded restriction
+  weights and coarse-level mean projection) and Dirichlet (inhomogeneous boundary data
+  supported). Grid dims must be 2^k+1 per active dimension. Grid-size-independent
+  convergence (< 0.15 residual reduction per V(2,2) cycle, 9²–129² validated); new
+  `POISSON_SOLVER_MG_SCALAR` convenience preset
+  (`lib/src/solvers/linear/cpu/linear_solver_multigrid.c`,
+  `lib/src/solvers/linear/cpu/multigrid_transfer.c`,
+  `tests/math/test_multigrid_operators.c`, `tests/math/test_multigrid_convergence.c`).
 - **Restart / checkpoint support** — portable, versioned, CRC-protected binary checkpoint
   format (`.cfdchk`) that saves and restores complete simulation state (grid, flow field,
   scalar params, time, solver name). Little-endian fixed-width encoding with an endianness
   marker and a format-version header that rejects unknown versions
   (`lib/src/io/checkpoint.c`, `lib/include/cfd/io/checkpoint.h`, `tests/io/test_checkpoint.c`).
 
+### Fixed
+
+- `poisson_solve_3d()` now checks `poisson_solver_init()`'s return status: a failed init
+  (e.g. multigrid on non-2^k+1 dims) no longer leaves a broken solver in the convenience
+  cache; the call returns -1 and later valid calls re-create the solver.
+
 ## [0.3.0] - 2026-06-23
 
 ### Added

CMakeLists.txt

@@ -351,6 +351,14 @@ if(BUILD_TESTS)
     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)
+    # mg_* transfer operators are not exported from the shared library, so the
+    # test must compile the translation unit directly there; in static builds
+    # they link from CFD::Library and compiling them again would duplicate the
+    # definitions.
+    add_executable(test_multigrid_operators
+        tests/math/test_multigrid_operators.c
+        $<$<BOOL:${BUILD_SHARED_LIBS}>:${CMAKE_CURRENT_SOURCE_DIR}/lib/src/solvers/linear/cpu/multigrid_transfer.c>)
+    add_executable(test_multigrid_convergence tests/math/test_multigrid_convergence.c)
 
     # Cross-architecture consistency tests
     add_executable(test_solver_architecture tests/validation/test_solver_architecture.c)
@@ -473,6 +481,8 @@ if(BUILD_TESTS)
     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>)
+    target_link_libraries(test_multigrid_operators PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
+    target_link_libraries(test_multigrid_convergence PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_solver_architecture PRIVATE CFD::Library unity)
 
     # Enable CTest support. include(CTest) generates DartConfiguration.tcl
@@ -635,6 +645,8 @@ if(BUILD_TESTS)
     add_test(NAME ResidualComputationTest COMMAND test_residual_computation)
     add_test(NAME CgScalingTest COMMAND test_cg_scaling)
     add_test(NAME NonuniformGridTest COMMAND test_nonuniform_grid)
+    add_test(NAME MultigridOperatorsTest COMMAND test_multigrid_operators)
+    add_test(NAME MultigridConvergenceTest COMMAND test_multigrid_convergence)
 
     # Label long-running validation tests to exclude from sanitizer job
     # These tests validate physics accuracy, not memory safety

ROADMAP.md

@@ -54,6 +54,7 @@ The single source of truth for backend gaps. Each algorithm targets scalar (CPU)
 |                     | CG / PCG       | done | done     | done     | done     | done |
 |                     | BiCGSTAB       | done | done     | done     | —        | done |
 |                     | GMRES(m)       | done | done     | done     | done     | —    |
+|                     | Multigrid (GMG)| done | —        | —        | —        | —    |
 | **Boundary Conds**  | All types      | done | done     | done     | done     | done |
 
 ### Known Limitations
@@ -73,6 +74,9 @@ Genuine constraints to be aware of (not backlog items):
 - **Modular library circular dependencies** — `cfd_scalar`/`cfd_simd` call `poisson_solve()`
   (in `cfd_api`) while `cfd_api` links against them. Resolved on Linux via linker groups;
   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.
 
 > The SOR optimal-ω and Jacobi spectral-radius formulas (ρ = cos(πh),
 > ω = 2/(1 + sin(πh))) assume Dirichlet BCs; with Neumann BCs optimal ω is typically lower
@@ -84,7 +88,7 @@ Genuine constraints to be aware of (not backlog items):
 
 | Phase | Theme | Priority | Status — what remains |
 | ----- | ----- | -------- | --------------------- |
-| 1 | Core Solver Improvements | P0–P3 | GMRES, multigrid, implicit integrators, SIMPLE/PISO, nonlinear & eigenvalue solvers |
+| 1 | Core Solver Improvements | P0–P3 | GMRES ✅, geometric multigrid ✅; AMG, implicit integrators, SIMPLE/PISO, nonlinear & eigenvalue solvers |
 | 2 | Physics Extensions | P1–P3 | RANS k-ε/SA ✅; realizable k-ε, k-ω SST remain; compressible, species, multiphase; energy-eq. extensions |
 | 3 | Geometry & Mesh | P1–P2 | Unstructured meshes, mesh I/O, adaptive refinement (3D ✅) |
 | 4 | Scalability & Performance | P1–P2 | MPI, GPU improvements, profiling tools (modular libs ✅) |
@@ -106,9 +110,9 @@ No-slip, Inlet, Outlet, Symmetry, Moving wall, Time-varying). See CHANGELOG.
 ### 1.2 Linear Solvers (P0)
 
 Implemented: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB, GMRES(m) (all with SIMD backends;
-CG is the default Poisson solver for projection methods). GPU standalone Jacobi, CG, Red-Black
-SOR, plain SOR, and BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses
-on-device CG.
+CG is the default Poisson solver for projection methods), and geometric multigrid (scalar).
+GPU standalone Jacobi, CG, Red-Black SOR, plain SOR, and BiCGSTAB are done and validated vs
+CPU; `solve_projection_method_gpu` uses on-device CG.
 
 GMRES(m) with restart is implemented across scalar, AVX2, NEON, and OMP backends (right-
 preconditioned Jacobi seam; validated vs CG on the SPD Poisson problem, with restart-no-stall
@@ -123,7 +127,11 @@ non-symmetric operators arrive, e.g. implicit advection-diffusion in §1.5).
 - [ ] GMRES GPU backend (deferred from the initial GMRES landing)
 - [ ] SSOR (Symmetric SOR) preconditioner
 - [ ] ILU preconditioner
-- [ ] Geometric multigrid
+- [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).
 - [ ] 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)
 

docs/guides/examples.md

@@ -684,7 +684,7 @@ CFD Platform Diagnostics
 1. Method comparison (Jacobi, SOR, Red-Black SOR, CG, CG+Jacobi PC, BiCGSTAB) on scalar backend
 2. Backend comparison (CG on Scalar, SIMD, OMP)
 3. Convenience API demo (`poisson_solve()`)
-4. Error handling (requesting unavailable multigrid solver)
+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.
 

docs/reference/api-reference.md

@@ -443,7 +443,7 @@ typedef enum {
     POISSON_METHOD_CG,            // Conjugate Gradient (SPD systems)
     POISSON_METHOD_BICGSTAB,      // BiCGSTAB (non-symmetric systems)
     POISSON_METHOD_GMRES,         // Restarted GMRES(m) (scalar/SIMD/OMP)
-    POISSON_METHOD_MULTIGRID      // Multigrid (future)
+    POISSON_METHOD_MULTIGRID      // Geometric multigrid V/W/F cycles (scalar; dims 2^k+1)
 } poisson_solver_method_t;
 ```
 
@@ -474,11 +474,24 @@ typedef struct {
     bool verbose;               // Print convergence info (default: false)
     poisson_precond_type_t preconditioner;  // Preconditioner (default: POISSON_PRECOND_NONE)
     int restart;                // GMRES(m) restart length (default: 0 = auto/30)
+
+    // Multigrid only (all 0-defaults are backward compatible)
+    mg_cycle_type_t mg_cycle;       // MG_CYCLE_V (default) / MG_CYCLE_W / MG_CYCLE_F
+    mg_smoother_type_t mg_smoother; // MG_SMOOTHER_REDBLACK_GS (default) / MG_SMOOTHER_JACOBI
+    mg_bc_type_t mg_bc;             // MG_BC_NEUMANN (default) / MG_BC_DIRICHLET
+    int mg_pre_smooth;              // Pre-smoothing sweeps (0 = default 2)
+    int mg_post_smooth;             // Post-smoothing sweeps (0 = default 2)
+    int mg_coarse_max_iter;         // Coarsest-grid smoother sweeps (0 = default 50)
+    int mg_max_levels;              // Max levels (0 = auto)
 } poisson_solver_params_t;
 
 poisson_solver_params_t poisson_solver_params_default(void);
 ```
 
+> Multigrid requires 2^k+1 grid points per active dimension (e.g. 33, 65, 129);
+> `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 Statistics
 
 ```c

docs/reference/solvers.md

@@ -255,6 +255,47 @@ poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_GMRES,
 poisson_solver_init(solver, nx, ny, nz, dx, dy, dz, &params);
 ```
 
+#### 8. Geometric Multigrid
+
+**Algorithm:** V-cycle (default), W-cycle, or F-cycle (full multigrid) on a
+hierarchy of coarsened grids. Each cycle pre-smooths, restricts the residual
+with full weighting, recursively solves the coarse correction equation,
+prolongates with bilinear (2D) / trilinear (3D) interpolation, and post-smooths.
+Smoothers: Red-Black Gauss-Seidel (default) or weighted Jacobi (ω=2/3).
+
+**Characteristics:**
+- **Grid-size-independent convergence**: ~0.05–0.1 residual reduction per
+  V(2,2) cycle regardless of resolution — O(N) total work, the optimal
+  complexity for structured-grid Poisson problems
+- **Grid constraint**: every active dimension must have 2^k+1 points
+  (5, 9, 17, 33, 65, 129, ...); other sizes return `CFD_ERROR_INVALID`
+- Two BC modes via `params.mg_bc`:
+  - `MG_BC_NEUMANN` (default) — zero-gradient BCs matching the other Poisson
+    solvers. The system is singular (solution defined up to a constant; RHS
+    should have zero interior mean). Restriction uses Neumann-folded boundary
+    weights and coarse RHS/corrections are mean-projected internally.
+  - `MG_BC_DIRICHLET` — caller-supplied boundary values of `x` are held fixed
+    (supports inhomogeneous data); coarse corrections use zero boundaries.
+- `MG_CYCLE_F` runs one full-multigrid pass (coarsest-first nested iteration)
+  on the first cycle — reaching discretization accuracy immediately — then
+  continues with V-cycles
+- 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.
+
+**Usage:**
+```c
+poisson_solver_params_t params = poisson_solver_params_default();
+params.mg_cycle = MG_CYCLE_V;        // or MG_CYCLE_W / MG_CYCLE_F
+params.mg_bc = MG_BC_NEUMANN;        // default; matches other solvers
+
+poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_MULTIGRID,
+                                                 POISSON_BACKEND_SCALAR);
+poisson_solver_init(solver, 65, 65, 1, dx, dy, 0.0, &params);  // dims 2^k+1
+```
+
 ### Linear Solver Performance Comparison
 
 **Problem:** 65×65 grid, tolerance = 1e-6
@@ -267,6 +308,7 @@ poisson_solver_init(solver, nx, ny, nz, dx, dy, dz, &params);
 | CG | ~80 | 5 | Best for large grids |
 | PCG (Jacobi) | ~80 | 5.5 | No benefit on uniform grid |
 | BiCGSTAB | ~40 | 4 | Fastest convergence |
+| Multigrid V(2,2) | ~8 cycles | — | O(N), grid-size-independent; needs 2^k+1 dims |
 
 **Note:** Jacobi preconditioning provides no benefit on uniform grids with constant coefficients (diagonal is constant 4/h²).
 

examples/poisson_solver_tuning.c

@@ -188,15 +188,17 @@ int main(void) {
         printf("  poisson_solve(CG_SCALAR): %d iterations, L2 error = %.2e\n", iters, err);
     }
 
-    /* Section 4: Error handling for unavailable solver */
+    /* Section 4: Error handling for unavailable solver.
+     * Multigrid only has a scalar backend, so requesting it on GPU always
+     * returns NULL (no silent fallbacks). */
     printf("\n--- Error Handling ---\n");
-    poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR);
+    poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_MULTIGRID, POISSON_BACKEND_GPU);
     if (!solver) {
-        printf("  Multigrid solver: not available (returns NULL)\n");
+        printf("  Multigrid GPU solver: not available (returns NULL)\n");
         printf("  Status: \"%s\"\n", cfd_get_error_string(cfd_get_last_status()));
         cfd_clear_error();
     } else {
-        printf("  Multigrid solver: available\n");
+        printf("  Multigrid GPU solver: available\n");
         poisson_solver_destroy(solver);
     }
 

lib/CMakeLists.txt

@@ -133,6 +133,8 @@ set(CFD_SCALAR_SOURCES
     src/solvers/linear/cpu/linear_solver_cg.c
     src/solvers/linear/cpu/linear_solver_bicgstab.c
     src/solvers/linear/cpu/linear_solver_gmres.c
+    src/solvers/linear/cpu/linear_solver_multigrid.c
+    src/solvers/linear/cpu/multigrid_transfer.c
 )
 
 # SIMD sources (AVX2/NEON optimized)

lib/include/cfd/solvers/poisson_solver.h

@@ -58,7 +58,7 @@ typedef enum {
     POISSON_METHOD_CG,            /**< Conjugate Gradient (for SPD systems) */
     POISSON_METHOD_BICGSTAB,      /**< BiCGSTAB (for non-symmetric systems) */
     POISSON_METHOD_GMRES,         /**< Restarted GMRES(m) (for non-symmetric systems) */
-    POISSON_METHOD_MULTIGRID      /**< Multigrid (future) */
+    POISSON_METHOD_MULTIGRID      /**< Geometric multigrid (V/W/F cycles, scalar backend) */
 } poisson_solver_method_t;
 
 /**
@@ -91,6 +91,38 @@ typedef enum {
     POISSON_PRECOND_JACOBI = 1  /**< Diagonal (Jacobi) preconditioning */
 } poisson_precond_type_t;
 
+/**
+ * Multigrid cycle type (POISSON_METHOD_MULTIGRID only)
+ */
+typedef enum {
+    MG_CYCLE_V = 0,   /**< V-cycle (default) */
+    MG_CYCLE_W = 1,   /**< W-cycle (two coarse-grid visits per level, more robust) */
+    MG_CYCLE_F = 2    /**< F-cycle / full multigrid (FMG first pass, then V-cycles) */
+} mg_cycle_type_t;
+
+/**
+ * Multigrid smoother type (POISSON_METHOD_MULTIGRID only)
+ */
+typedef enum {
+    MG_SMOOTHER_REDBLACK_GS = 0,  /**< Red-Black Gauss-Seidel, omega=1 (default) */
+    MG_SMOOTHER_JACOBI = 1        /**< Weighted Jacobi, omega=2/3 */
+} mg_smoother_type_t;
+
+/**
+ * Multigrid boundary-condition mode (POISSON_METHOD_MULTIGRID only)
+ *
+ * MG_BC_NEUMANN matches the default zero-gradient behavior of all other
+ * Poisson solvers (solution defined up to an additive constant; the RHS
+ * should have zero interior mean for full convergence).
+ * MG_BC_DIRICHLET holds the caller-supplied boundary values of x fixed on
+ * the finest grid and uses homogeneous zero boundaries for coarse-level
+ * corrections (supports inhomogeneous Dirichlet data).
+ */
+typedef enum {
+    MG_BC_NEUMANN = 0,   /**< Zero-gradient BCs (default; matches other solvers) */
+    MG_BC_DIRICHLET = 1  /**< Fixed boundary values on finest level */
+} mg_bc_type_t;
+
 /* ============================================================================
  * PARAMETERS AND STATISTICS
  * ============================================================================ */
@@ -107,6 +139,15 @@ typedef struct {
     bool verbose;              /**< Print iteration progress (default: false) */
     poisson_precond_type_t preconditioner; /**< Preconditioner type (default: NONE) */
     int restart;               /**< GMRES restart length m (default: 0 = auto/30); ignored by other methods */
+
+    /* Multigrid parameters (POISSON_METHOD_MULTIGRID only; 0 = backward-compatible default) */
+    mg_cycle_type_t mg_cycle;       /**< Cycle type (default: MG_CYCLE_V) */
+    mg_smoother_type_t mg_smoother; /**< Smoother (default: MG_SMOOTHER_REDBLACK_GS) */
+    mg_bc_type_t mg_bc;             /**< Boundary-condition mode (default: MG_BC_NEUMANN) */
+    int mg_pre_smooth;              /**< Pre-smoothing sweeps per level (0 = default 2) */
+    int mg_post_smooth;             /**< Post-smoothing sweeps per level (0 = default 2) */
+    int mg_coarse_max_iter;         /**< Smoother sweeps on coarsest grid (0 = default 50) */
+    int mg_max_levels;              /**< Max grid levels (0 = auto: coarsen as far as possible) */
 } poisson_solver_params_t;
 
 /**
@@ -130,6 +171,8 @@ typedef struct {
  * - omega: 0.0 (auto-compute optimal for grid dimensions)
  * - check_interval: 1
  * - verbose: false
+ * - mg_cycle: MG_CYCLE_V, mg_smoother: MG_SMOOTHER_REDBLACK_GS, mg_bc: MG_BC_NEUMANN
+ * - mg_pre_smooth/mg_post_smooth: 0 (= 2), mg_coarse_max_iter: 0 (= 50), mg_max_levels: 0 (= auto)
  */
 CFD_LIBRARY_EXPORT poisson_solver_params_t poisson_solver_params_default(void);
 
@@ -396,6 +439,7 @@ CFD_LIBRARY_EXPORT bool poisson_solver_backend_available(poisson_solver_backend_
 #define POISSON_SOLVER_TYPE_GMRES_SCALAR      "gmres_scalar"
 #define POISSON_SOLVER_TYPE_GMRES_OMP         "gmres_omp"
 #define POISSON_SOLVER_TYPE_GMRES_SIMD        "gmres_simd"
+#define POISSON_SOLVER_TYPE_MG_SCALAR         "multigrid_scalar"
 
 /* ============================================================================
  * CONVENIENCE API
@@ -416,7 +460,8 @@ typedef enum {
     POISSON_SOLVER_CG_SCALAR = 5,      /**< Conjugate Gradient with scalar backend (always available) */
     POISSON_SOLVER_CG_SIMD = 6,        /**< Conjugate Gradient with SIMD backend (runtime detection) */
     POISSON_SOLVER_CG_OMP = 7,         /**< Conjugate Gradient with OpenMP backend */
-    POISSON_SOLVER_SOR_SIMD = 8        /**< SOR with SIMD backend (Block SOR, runtime detection) */
+    POISSON_SOLVER_SOR_SIMD = 8,       /**< SOR with SIMD backend (Block SOR, runtime detection) */
+    POISSON_SOLVER_MG_SCALAR = 9       /**< Geometric multigrid with scalar backend (grid dims must be 2^k+1) */
 } poisson_solver_type;
 
 /** Default Poisson solver - uses runtime SIMD detection */