From 15b04a992d096109cbd6367bd031635bd5fe9ec7 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 11 Jul 2026 20:50:57 +0300 Subject: [PATCH 1/4] Add geometric multigrid Poisson solver (scalar backend) CG iteration counts grow O(N) with grid refinement; geometric multigrid converges in a grid-size-independent number of cycles for O(N) total work, making it the optimal solver for large structured-grid pressure Poisson problems (ROADMAP 1.2). V/W/F(FMG) cycles with Red-Black Gauss-Seidel or weighted-Jacobi smoothers, full-weighting restriction and bilinear/trilinear prolongation, 2D/3D. Grid dims must be 2^k+1 per active dimension. Two BC modes, because the spec's Dirichlet assumption did not match the codebase: MG_BC_NEUMANN (default) mirrors the zero-gradient convention of the other solvers, MG_BC_DIRICHLET holds caller boundary values fixed. The Neumann mirror slaves the boundary vertex to its interior neighbor, so restriction folds the toward-boundary weights to stay the exact adjoint of prolongation - without the fold the convergence factor degrades several-fold. The singular Neumann system is handled by mean-projecting restricted RHS vectors and coarse corrections; the coarsest Neumann grid stops at 5 points per dim (a 3x3 mirror grid has an identically zero operator). Neumann cross-checks use RB-SOR, not CG: CG freezes boundary values during Krylov iterations (BC applied only at start/end) and therefore solves a different discrete system than the per-sweep-BC solvers. Also fix poisson_solve_3d ignoring init status - a failed init (now possible via the 2^k+1 check) no longer caches a broken solver. Verified: V(2,2) factor < 0.15, 9x9-129x129 converge in <= 15 cycles (spread <= 3), FMG reaches discretization accuracy in one pass; full suite 103/103 incl. Ghia/cavity validation. --- CHANGELOG.md | 17 + CMakeLists.txt | 8 + ROADMAP.md | 18 +- docs/guides/examples.md | 2 +- docs/reference/api-reference.md | 15 +- docs/reference/solvers.md | 42 + examples/poisson_solver_tuning.c | 10 +- lib/CMakeLists.txt | 2 + lib/include/cfd/solvers/poisson_solver.h | 49 +- .../linear/cpu/linear_solver_multigrid.c | 575 +++++++++++++ .../solvers/linear/cpu/multigrid_transfer.c | 200 +++++ lib/src/solvers/linear/linear_solver.c | 40 +- .../solvers/linear/linear_solver_internal.h | 3 + lib/src/solvers/linear/multigrid_internal.h | 125 +++ tests/math/test_multigrid_convergence.c | 783 ++++++++++++++++++ tests/math/test_multigrid_operators.c | 681 +++++++++++++++ 16 files changed, 2553 insertions(+), 17 deletions(-) create mode 100644 lib/src/solvers/linear/cpu/linear_solver_multigrid.c create mode 100644 lib/src/solvers/linear/cpu/multigrid_transfer.c create mode 100644 lib/src/solvers/linear/multigrid_internal.h create mode 100644 tests/math/test_multigrid_convergence.c create mode 100644 tests/math/test_multigrid_operators.c diff --git a/CHANGELOG.md b/CHANGELOG.md index 04443853..74eca05c 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 97f68010..ac04e8c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -351,6 +351,10 @@ 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) + add_executable(test_multigrid_operators + tests/math/test_multigrid_operators.c + 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 +477,8 @@ if(BUILD_TESTS) target_link_libraries(test_residual_computation PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_cg_scaling PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_nonuniform_grid PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_multigrid_operators PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_multigrid_convergence PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_solver_architecture PRIVATE CFD::Library unity) # Enable CTest support. include(CTest) generates DartConfiguration.tcl @@ -635,6 +641,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 diff --git a/ROADMAP.md b/ROADMAP.md index 8a5242f4..587d1a50 100644 --- a/ROADMAP.md +++ b/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) diff --git a/docs/guides/examples.md b/docs/guides/examples.md index 51032d1e..413b6491 100644 --- a/docs/guides/examples.md +++ b/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. diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md index a444a125..b25484ed 100644 --- a/docs/reference/api-reference.md +++ b/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 diff --git a/docs/reference/solvers.md b/docs/reference/solvers.md index c7ebf083..037a944b 100644 --- a/docs/reference/solvers.md +++ b/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, ¶ms); ``` +#### 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, ¶ms); // 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, ¶ms); | 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²). diff --git a/examples/poisson_solver_tuning.c b/examples/poisson_solver_tuning.c index 7c1749aa..5d507871 100644 --- a/examples/poisson_solver_tuning.c +++ b/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); } diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 665d2ddb..b0ea9557 100644 --- a/lib/CMakeLists.txt +++ b/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) diff --git a/lib/include/cfd/solvers/poisson_solver.h b/lib/include/cfd/solvers/poisson_solver.h index b7887969..40d892e2 100644 --- a/lib/include/cfd/solvers/poisson_solver.h +++ b/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 */ diff --git a/lib/src/solvers/linear/cpu/linear_solver_multigrid.c b/lib/src/solvers/linear/cpu/linear_solver_multigrid.c new file mode 100644 index 00000000..5b63c097 --- /dev/null +++ b/lib/src/solvers/linear/cpu/linear_solver_multigrid.c @@ -0,0 +1,575 @@ +/** + * @file linear_solver_multigrid.c + * @brief Geometric multigrid Poisson solver - scalar CPU implementation + * + * V/W/F(FMG) cycles with Red-Black Gauss-Seidel (default) or weighted Jacobi + * smoothing, full-weighting restriction and bilinear/trilinear prolongation + * (see cpu/multigrid_transfer.c). Grid dimensions must be 2^k+1 per active + * dimension. + * + * Boundary-condition modes (params.mg_bc): + * - MG_BC_NEUMANN (default): zero-gradient BCs matching the other Poisson + * solvers. The system is singular (constant nullspace); restricted RHS + * vectors are projected to zero interior mean, coarse corrections are + * mean-subtracted, and the solution converges up to an additive constant. + * - MG_BC_DIRICHLET: caller-supplied boundary values of x are held fixed on + * the finest grid; coarse-level corrections use homogeneous zero boundaries. + * + * Level 0 (finest) borrows the caller's x/rhs arrays; coarser levels own + * their buffers. The public x_temp argument is unused (like CG); the Jacobi + * smoother allocates its own per-level temporaries. + */ + +#include "../linear_solver_internal.h" +#include "../multigrid_internal.h" + +#include "cfd/boundary/boundary_conditions.h" +#include "cfd/core/indexing.h" +#include "cfd/core/memory.h" + +#include + +/* ============================================================================ + * BOUNDARY CONDITIONS + * ============================================================================ */ + +/** + * Apply the BC mode to a field at one level. + * + * Neumann: mirror of poisson_solver_apply_bc()'s default (z-plane copies plus + * per-plane zero-gradient), but with this level's dimensions. + * Dirichlet: no-op — boundary values are never touched anywhere in MG, so + * finest-level data stays fixed and coarse corrections stay zero. + */ +static void mg_apply_bc_level(const mg_level_t* L, mg_bc_type_t bc, double* x) { + if (bc == MG_BC_DIRICHLET) { + return; + } + + size_t nx = L->nx; + size_t ny = L->ny; + size_t nz = L->nz; + size_t plane_size = nx * ny; + + if (nz > 1) { + memcpy(x, x + plane_size, plane_size * sizeof(double)); + memcpy(x + (nz - 1) * plane_size, + x + (nz - 2) * plane_size, + plane_size * sizeof(double)); + } + + for (size_t k = 0; k < nz; k++) { + bc_apply_scalar_cpu(x + k * plane_size, nx, ny, BC_TYPE_NEUMANN); + } +} + +/* ============================================================================ + * SMOOTHERS + * ============================================================================ */ + +/** + * Red-Black Gauss-Seidel smoother (omega = 1: over-relaxation degrades the + * high-frequency smoothing factor, so plain GS is correct inside multigrid). + * Same stencil algebra as linear_solver_redblack.c. BCs after each sweep. + */ +static void mg_smooth_rbgs(const mg_level_t* L, mg_bc_type_t bc, + double* x, const double* rhs, int sweeps) { + size_t nx = L->nx; + size_t ny = L->ny; + double dx2 = L->dx2; + double dy2 = L->dy2; + double inv_dz2 = L->inv_dz2; + double inv_factor = L->inv_factor; + size_t stride_z = L->stride_z; + + for (int sweep = 0; sweep < sweeps; sweep++) { + /* Red pass: (i+j+k) % 2 == 0, then black pass: (i+j+k) % 2 == 1 */ + for (int color = 0; color < 2; color++) { + for (size_t k = L->k_start; k < L->k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + size_t i_start = + ((j + k) % 2 == (size_t)color) ? 1 : 2; + for (size_t i = i_start; i < nx - 1; i += 2) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + + x[idx] = -(rhs[idx] + - (x[idx + 1] + x[idx - 1]) / dx2 + - (x[idx + nx] + x[idx - nx]) / dy2 + - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2 + ) * inv_factor; + } + } + } + } + + mg_apply_bc_level(L, bc, x); + } +} + +/** + * Weighted Jacobi smoother (omega = 2/3: optimal high-frequency damping). + * Copies interior points only from x_temp back to x, preserving boundary + * values in Dirichlet mode. BCs after each sweep. + */ +static void mg_smooth_jacobi(const mg_level_t* L, mg_bc_type_t bc, + double* x, double* x_temp, const double* rhs, + int sweeps) { + const double omega = 2.0 / 3.0; + size_t nx = L->nx; + size_t ny = L->ny; + double dx2 = L->dx2; + double dy2 = L->dy2; + double inv_dz2 = L->inv_dz2; + double inv_factor = L->inv_factor; + size_t stride_z = L->stride_z; + + for (int sweep = 0; sweep < sweeps; sweep++) { + for (size_t k = L->k_start; k < L->k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + + double x_jacobi = -(rhs[idx] + - (x[idx + 1] + x[idx - 1]) / dx2 + - (x[idx + nx] + x[idx - nx]) / dy2 + - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2 + ) * inv_factor; + + x_temp[idx] = x[idx] + omega * (x_jacobi - x[idx]); + } + } + } + + for (size_t k = L->k_start; k < L->k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + x[idx] = x_temp[idx]; + } + } + } + + mg_apply_bc_level(L, bc, x); + } +} + +static void mg_smooth(const mg_context_t* ctx, const mg_level_t* L, + double* x, const double* rhs, int sweeps) { + if (ctx->smoother_type == MG_SMOOTHER_JACOBI) { + mg_smooth_jacobi(L, ctx->bc_mode, x, L->x_temp, rhs, sweeps); + } else { + mg_smooth_rbgs(L, ctx->bc_mode, x, rhs, sweeps); + } +} + +/* ============================================================================ + * RESIDUAL AND GRID TRANSFERS + * ============================================================================ */ + +/** + * r = rhs - Laplacian(x) at interior points. The coarse problem is + * Laplacian(e) = r, so that Laplacian(x + P e) ~ rhs. Boundary entries of r + * are never written (they stay zero from allocation; restriction never + * reads them). + */ +static void mg_residual(const mg_level_t* L, const double* x, + const double* rhs, double* r) { + size_t nx = L->nx; + size_t ny = L->ny; + double dx2 = L->dx2; + double dy2 = L->dy2; + double inv_dz2 = L->inv_dz2; + size_t stride_z = L->stride_z; + + for (size_t k = L->k_start; k < L->k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D(i, j, nx); + + double laplacian = + (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) / dx2 + + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) / dy2 + + (x[idx + stride_z] + x[idx - stride_z] + - 2.0 * x[idx]) * inv_dz2; + + r[idx] = rhs[idx] - laplacian; + } + } + } +} + +/** + * Restrict a fine-level field into a coarse-level (MG-owned) RHS buffer. + * In Neumann mode the coarse RHS is projected to zero interior mean — the + * singular mirror-BC operator only admits solutions for compatible RHS. + */ +static void mg_restrict_level(const mg_context_t* ctx, + const mg_level_t* Lf, const mg_level_t* Lc, + const double* fine, double* coarse) { + int fold = (ctx->bc_mode == MG_BC_NEUMANN); + + if (Lf->nz > 1) { + mg_restrict_3d(fine, coarse, Lf->nx, Lf->ny, Lf->nz, + Lc->nx, Lc->ny, Lc->nz, fold); + } else { + mg_restrict_2d(fine, coarse, Lf->nx, Lf->ny, Lc->nx, Lc->ny, fold); + } + + if (ctx->bc_mode == MG_BC_NEUMANN) { + mg_subtract_interior_mean(coarse, Lc->nx, Lc->ny, Lc->nz); + } +} + +static void mg_prolongate_level(const mg_level_t* Lc, const mg_level_t* Lf, + const double* coarse, double* fine) { + if (Lf->nz > 1) { + mg_prolongate_add_3d(coarse, fine, Lc->nx, Lc->ny, Lc->nz, + Lf->nx, Lf->ny, Lf->nz); + } else { + mg_prolongate_add_2d(coarse, fine, Lc->nx, Lc->ny, Lf->nx, Lf->ny); + } +} + +/* ============================================================================ + * CYCLES + * ============================================================================ */ + +/** + * Coarsest-grid solve: over-smoothing. A 3x3 Dirichlet grid (one unknown) is + * solved exactly by the first sweep; Neumann coarsest grids (>= 5x5, RHS + * already mean-projected) converge geometrically in a few dozen sweeps at + * negligible cost. + */ +static void mg_coarse_solve(const mg_context_t* ctx, const mg_level_t* L, + double* x, const double* rhs) { + mg_smooth(ctx, L, x, rhs, ctx->coarse_max_iter); +} + +/** Recursive V/W cycle. Level 0 is called with the caller's arrays. */ +static void mg_cycle(mg_context_t* ctx, int level, double* x, + const double* rhs) { + mg_level_t* L = &ctx->levels[level]; + + if (level == ctx->num_levels - 1) { + mg_coarse_solve(ctx, L, x, rhs); + return; + } + + mg_level_t* Lc = &ctx->levels[level + 1]; + + /* Pre-smooth, then form and restrict the residual */ + mg_smooth(ctx, L, x, rhs, ctx->nu1); + mg_residual(L, x, rhs, L->residual); + mg_restrict_level(ctx, L, Lc, L->residual, Lc->rhs); + + /* Coarse-grid correction with zero initial guess. W-cycle: the second + * visit continues refining the same correction. */ + memset(Lc->x, 0, Lc->total * sizeof(double)); + int num_visits = (ctx->cycle_type == MG_CYCLE_W) ? 2 : 1; + for (int v = 0; v < num_visits; v++) { + mg_cycle(ctx, level + 1, Lc->x, Lc->rhs); + } + + if (ctx->bc_mode == MG_BC_NEUMANN) { + /* Constants are in the operator's nullspace: removing the mean keeps + * the fine solution from drifting. Mirror the boundary so + * prolongation interpolates consistent ghost values. */ + mg_subtract_interior_mean(Lc->x, Lc->nx, Lc->ny, Lc->nz); + mg_apply_bc_level(Lc, ctx->bc_mode, Lc->x); + } + + mg_prolongate_level(Lc, L, Lc->x, x); + mg_apply_bc_level(L, ctx->bc_mode, x); + + mg_smooth(ctx, L, x, rhs, ctx->nu2); +} + +/** + * FMG (full multigrid) nested iteration: restrict the original RHS through + * all levels, solve on the coarsest, then work upward prolongating the + * solution as the initial guess and running one V-cycle per level. Replaces + * the caller's initial guess (interior only; boundary data preserved). + * + * Ascent ordering is load-bearing: levels[l+1].x must be consumed before + * mg_cycle(l) overwrites levels[l+1].{x,rhs} with residual-correction data. + */ +static void mg_fmg(mg_context_t* ctx, double* x, const double* rhs) { + int nl = ctx->num_levels; + mg_level_t* levels = ctx->levels; + + if (nl == 1) { + mg_cycle(ctx, 0, x, rhs); + return; + } + + /* Descent: levels[l].rhs = restricted original b for l >= 1 */ + mg_restrict_level(ctx, &levels[0], &levels[1], rhs, levels[1].rhs); + for (int l = 1; l <= nl - 2; l++) { + mg_restrict_level(ctx, &levels[l], &levels[l + 1], + levels[l].rhs, levels[l + 1].rhs); + } + + /* Coarsest solve */ + memset(levels[nl - 1].x, 0, levels[nl - 1].total * sizeof(double)); + mg_coarse_solve(ctx, &levels[nl - 1], levels[nl - 1].x, + levels[nl - 1].rhs); + + /* Ascent through MG-owned levels */ + for (int l = nl - 2; l >= 1; l--) { + mg_apply_bc_level(&levels[l + 1], ctx->bc_mode, levels[l + 1].x); + memset(levels[l].x, 0, levels[l].total * sizeof(double)); + mg_prolongate_level(&levels[l + 1], &levels[l], + levels[l + 1].x, levels[l].x); + mg_apply_bc_level(&levels[l], ctx->bc_mode, levels[l].x); + mg_cycle(ctx, l, levels[l].x, levels[l].rhs); + } + + /* Finest level: zero the interior of the caller's x (boundary data must + * survive in Dirichlet mode), prolongate the level-1 solution in, and + * finish with one V-cycle */ + mg_apply_bc_level(&levels[1], ctx->bc_mode, levels[1].x); + { + mg_level_t* L0 = &levels[0]; + for (size_t k = L0->k_start; k < L0->k_end; k++) { + for (size_t j = 1; j < L0->ny - 1; j++) { + for (size_t i = 1; i < L0->nx - 1; i++) { + x[k * L0->stride_z + IDX_2D(i, j, L0->nx)] = 0.0; + } + } + } + } + mg_prolongate_level(&levels[1], &levels[0], levels[1].x, x); + mg_apply_bc_level(&levels[0], ctx->bc_mode, x); + mg_cycle(ctx, 0, x, rhs); +} + +/* ============================================================================ + * LIFECYCLE + * ============================================================================ */ + +static void mg_free_levels(mg_context_t* ctx) { + if (!ctx->levels) { + return; + } + for (int l = 0; l < ctx->num_levels; l++) { + cfd_free(ctx->levels[l].x); + cfd_free(ctx->levels[l].x_temp); + cfd_free(ctx->levels[l].rhs); + cfd_free(ctx->levels[l].residual); + } + cfd_free(ctx->levels); + ctx->levels = NULL; + ctx->num_levels = 0; +} + +static void mg_scalar_destroy(poisson_solver_t* solver) { + if (solver && solver->context) { + mg_context_t* ctx = (mg_context_t*)solver->context; + mg_free_levels(ctx); + cfd_free(ctx); + solver->context = NULL; + } +} + +static cfd_status_t mg_scalar_init( + poisson_solver_t* solver, + size_t nx, size_t ny, size_t nz, + double dx, double dy, double dz, + const poisson_solver_params_t* params) +{ + /* Geometric coarsening requires 2^k+1 points per active dimension */ + if (!mg_is_pow2_plus1(nx) || !mg_is_pow2_plus1(ny) || + (nz > 1 && !mg_is_pow2_plus1(nz))) { + return CFD_ERROR_INVALID; + } + + if ((int)params->mg_cycle < MG_CYCLE_V || + (int)params->mg_cycle > MG_CYCLE_F || + (int)params->mg_smoother < MG_SMOOTHER_REDBLACK_GS || + (int)params->mg_smoother > MG_SMOOTHER_JACOBI || + (int)params->mg_bc < MG_BC_NEUMANN || + (int)params->mg_bc > MG_BC_DIRICHLET || + params->mg_pre_smooth < 0 || params->mg_post_smooth < 0 || + params->mg_coarse_max_iter < 0 || params->mg_max_levels < 0) { + return CFD_ERROR_INVALID; + } + + /* Re-init support: drop any previous hierarchy */ + mg_scalar_destroy(solver); + + mg_context_t* ctx = (mg_context_t*)cfd_calloc(1, sizeof(mg_context_t)); + if (!ctx) { + return CFD_ERROR_NOMEM; + } + + ctx->cycle_type = params->mg_cycle; + ctx->smoother_type = params->mg_smoother; + ctx->bc_mode = params->mg_bc; + ctx->nu1 = (params->mg_pre_smooth > 0) + ? params->mg_pre_smooth : MG_DEFAULT_PRE_SMOOTH; + ctx->nu2 = (params->mg_post_smooth > 0) + ? params->mg_post_smooth : MG_DEFAULT_POST_SMOOTH; + ctx->coarse_max_iter = (params->mg_coarse_max_iter > 0) + ? params->mg_coarse_max_iter : MG_DEFAULT_COARSE_MAX_ITER; + + /* Count levels: coarsen all active dimensions simultaneously while every + * next dimension stays at or above the BC-mode floor */ + size_t min_coarse = (ctx->bc_mode == MG_BC_NEUMANN) + ? MG_MIN_COARSE_DIM_NEUMANN : MG_MIN_COARSE_DIM_DIRICHLET; + int num_levels = 1; + { + size_t cx = nx, cy = ny, cz = nz; + while (params->mg_max_levels == 0 || + num_levels < params->mg_max_levels) { + size_t nx2 = (cx - 1) / 2 + 1; + size_t ny2 = (cy - 1) / 2 + 1; + size_t nz2 = (nz > 1) ? (cz - 1) / 2 + 1 : 1; + if (nx2 < min_coarse || ny2 < min_coarse || + (nz > 1 && nz2 < min_coarse)) { + break; + } + cx = nx2; + cy = ny2; + cz = nz2; + num_levels++; + } + } + ctx->num_levels = num_levels; + + ctx->levels = (mg_level_t*)cfd_calloc((size_t)num_levels, + sizeof(mg_level_t)); + if (!ctx->levels) { + cfd_free(ctx); + return CFD_ERROR_NOMEM; + } + + for (int l = 0; l < num_levels; l++) { + mg_level_t* L = &ctx->levels[l]; + double scale = (double)((size_t)1 << l); + + L->nx = ((nx - 1) >> l) + 1; + L->ny = ((ny - 1) >> l) + 1; + L->nz = (nz > 1) ? ((nz - 1) >> l) + 1 : 1; + L->total = L->nx * L->ny * L->nz; + L->dx2 = (dx * scale) * (dx * scale); + L->dy2 = (dy * scale) * (dy * scale); + L->inv_dz2 = poisson_solver_compute_inv_dz2(dz * scale); + poisson_solver_compute_3d_bounds(L->nz, L->nx, L->ny, + &L->stride_z, &L->k_start, &L->k_end); + L->inv_factor = 1.0 / (2.0 / L->dx2 + 2.0 / L->dy2 + 2.0 * L->inv_dz2); + + /* Level 0 borrows the caller's x/rhs; the coarsest level needs no + * residual buffer */ + int need_owned = (l >= 1); + int need_residual = (l <= num_levels - 2); + int need_temp = (ctx->smoother_type == MG_SMOOTHER_JACOBI); + + if ((need_owned && + (!(L->x = (double*)cfd_calloc(L->total, sizeof(double))) || + !(L->rhs = (double*)cfd_calloc(L->total, sizeof(double))))) || + (need_residual && + !(L->residual = (double*)cfd_calloc(L->total, sizeof(double)))) || + (need_temp && + !(L->x_temp = (double*)cfd_calloc(L->total, sizeof(double))))) { + mg_free_levels(ctx); + cfd_free(ctx); + return CFD_ERROR_NOMEM; + } + } + + ctx->fmg_pending = (ctx->cycle_type == MG_CYCLE_F); + ctx->initialized = 1; + solver->context = ctx; + return CFD_SUCCESS; +} + +/* ============================================================================ + * SOLVE / ITERATE / APPLY_BC + * ============================================================================ */ + +static cfd_status_t mg_scalar_iterate( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + double* residual) +{ + (void)x_temp; /* MG allocates its own temporaries */ + + mg_context_t* ctx = (mg_context_t*)solver->context; + if (!ctx || !ctx->initialized) { + return CFD_ERROR_INVALID; + } + + if (ctx->fmg_pending) { + mg_fmg(ctx, x, rhs); + ctx->fmg_pending = 0; + } else { + mg_cycle(ctx, 0, x, rhs); + } + + if (residual) { + *residual = poisson_solver_compute_residual(solver, x, rhs); + } + + return CFD_SUCCESS; +} + +/** + * Custom solve wrapper: the FMG flag must reset at the start of every solve, + * and iterate() has no solve-start hook. Delegates to the common loop (each + * "iteration" is one cycle). + */ +static cfd_status_t mg_scalar_solve( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + poisson_solver_stats_t* stats) +{ + mg_context_t* ctx = (mg_context_t*)solver->context; + if (!ctx || !ctx->initialized) { + return CFD_ERROR_INVALID; + } + + ctx->fmg_pending = (ctx->cycle_type == MG_CYCLE_F); + return poisson_solver_solve_common(solver, x, x_temp, rhs, stats); +} + +/** + * Non-NULL apply_bc is required: the public poisson_solver_apply_bc() default + * is Neumann, which would corrupt Dirichlet boundary data. + */ +static void mg_scalar_apply_bc(poisson_solver_t* solver, double* x) { + mg_context_t* ctx = (mg_context_t*)solver->context; + if (!ctx || !ctx->initialized) { + return; + } + mg_apply_bc_level(&ctx->levels[0], ctx->bc_mode, x); +} + +/* ============================================================================ + * FACTORY FUNCTION + * ============================================================================ */ + +poisson_solver_t* create_multigrid_scalar_solver(void) { + poisson_solver_t* solver = + (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t)); + if (!solver) { + return NULL; + } + + solver->name = POISSON_SOLVER_TYPE_MG_SCALAR; + solver->description = "Geometric multigrid V/W/F-cycle (scalar CPU)"; + solver->method = POISSON_METHOD_MULTIGRID; + solver->backend = POISSON_BACKEND_SCALAR; + solver->params = poisson_solver_params_default(); + + solver->init = mg_scalar_init; + solver->destroy = mg_scalar_destroy; + solver->solve = mg_scalar_solve; /* Resets FMG state, then common loop */ + solver->iterate = mg_scalar_iterate; /* One V/W cycle (or the FMG pass) */ + solver->apply_bc = mg_scalar_apply_bc; + + return solver; +} diff --git a/lib/src/solvers/linear/cpu/multigrid_transfer.c b/lib/src/solvers/linear/cpu/multigrid_transfer.c new file mode 100644 index 00000000..ce351dde --- /dev/null +++ b/lib/src/solvers/linear/cpu/multigrid_transfer.c @@ -0,0 +1,200 @@ +/** + * @file multigrid_transfer.c + * @brief Grid-transfer operators for geometric multigrid (restriction/prolongation) + * + * Pure array math with no library runtime dependencies, so the operator unit + * test can compile this translation unit directly (see multigrid_internal.h). + * + * Grid relationship: nxf = 2*nxc - 1 (fine dims are 2^k+1, coarse (2^(k-1))+1). + * Coarse interior point (I,J,K) sits at fine point (2I,2J,2K); its + * full-weighting stencil spans fine indices 2I-1..2I+1, which always lie in + * the fine interior for coarse-interior I. Restriction therefore never reads + * fine boundary values, and both operators write interior points only. + */ + +#include "../multigrid_internal.h" + +#include "cfd/core/indexing.h" + +/* ============================================================================ + * RESTRICTION (fine -> coarse, full weighting) + * ============================================================================ */ + +/** + * Per-dimension 1D restriction weights at a coarse index: (1,2,1) in the + * interior. With fold_neumann, the toward-boundary weight doubles at + * boundary-adjacent coarse points (I==1 and/or I==n_coarse-2) — the exact + * adjoint of prolongation through mirrored (zero-gradient) ghost values. + */ +static void mg_weights_1d(size_t idx, size_t n_coarse, int fold_neumann, + double w[3]) { + w[0] = (fold_neumann && idx == 1) ? 2.0 : 1.0; + w[1] = 2.0; + w[2] = (fold_neumann && idx == n_coarse - 2) ? 2.0 : 1.0; +} + +void mg_restrict_2d(const double* fine, double* coarse, + size_t nxf, size_t nyf, size_t nxc, size_t nyc, + int fold_neumann) { + (void)nyf; + /* Separable full weighting: interior stencil (1/16)*(4 center, 2 edge, + * 1 corner) */ + for (size_t J = 1; J < nyc - 1; J++) { + double wy[3]; + mg_weights_1d(J, nyc, fold_neumann, wy); + for (size_t I = 1; I < nxc - 1; I++) { + double wx[3]; + mg_weights_1d(I, nxc, fold_neumann, wx); + + size_t i = 2 * I; + size_t j = 2 * J; + double sum = 0.0; + for (int oj = -1; oj <= 1; oj++) { + for (int oi = -1; oi <= 1; oi++) { + double w = wx[oi + 1] * wy[oj + 1]; + sum += w * fine[IDX_2D(i + oi, j + oj, nxf)]; + } + } + coarse[IDX_2D(I, J, nxc)] = sum / 16.0; + } + } +} + +void mg_restrict_3d(const double* fine, double* coarse, + size_t nxf, size_t nyf, size_t nzf, + size_t nxc, size_t nyc, size_t nzc, + int fold_neumann) { + (void)nzf; + /* Separable 27-point full weighting: interior stencil (1/64)*(8 center, + * 4 face, 2 edge, 1 corner) */ + for (size_t K = 1; K < nzc - 1; K++) { + double wz[3]; + mg_weights_1d(K, nzc, fold_neumann, wz); + for (size_t J = 1; J < nyc - 1; J++) { + double wy[3]; + mg_weights_1d(J, nyc, fold_neumann, wy); + for (size_t I = 1; I < nxc - 1; I++) { + double wx[3]; + mg_weights_1d(I, nxc, fold_neumann, wx); + + size_t i = 2 * I; + size_t j = 2 * J; + size_t k = 2 * K; + double sum = 0.0; + for (int ok = -1; ok <= 1; ok++) { + for (int oj = -1; oj <= 1; oj++) { + for (int oi = -1; oi <= 1; oi++) { + double w = wx[oi + 1] * wy[oj + 1] * wz[ok + 1]; + sum += w * fine[IDX_3D(i + oi, j + oj, k + ok, + nxf, nyf)]; + } + } + } + coarse[IDX_3D(I, J, K, nxc, nyc)] = sum / 64.0; + } + } + } +} + +/* ============================================================================ + * PROLONGATION (coarse -> fine, bilinear/trilinear, additive) + * ============================================================================ */ + +void mg_prolongate_add_2d(const double* coarse, double* fine, + size_t nxc, size_t nyc, size_t nxf, size_t nyf) { + (void)nyc; + for (size_t j = 1; j < nyf - 1; j++) { + size_t J = j >> 1; + int jodd = (int)(j & 1); + for (size_t i = 1; i < nxf - 1; i++) { + size_t I = i >> 1; + int iodd = (int)(i & 1); + double e; + if (!iodd && !jodd) { + e = coarse[IDX_2D(I, J, nxc)]; + } else if (iodd && !jodd) { + e = 0.5 * (coarse[IDX_2D(I, J, nxc)] + + coarse[IDX_2D(I + 1, J, nxc)]); + } else if (!iodd && jodd) { + e = 0.5 * (coarse[IDX_2D(I, J, nxc)] + + coarse[IDX_2D(I, J + 1, nxc)]); + } else { + e = 0.25 * (coarse[IDX_2D(I, J, nxc)] + + coarse[IDX_2D(I + 1, J, nxc)] + + coarse[IDX_2D(I, J + 1, nxc)] + + coarse[IDX_2D(I + 1, J + 1, nxc)]); + } + fine[IDX_2D(i, j, nxf)] += e; + } + } +} + +void mg_prolongate_add_3d(const double* coarse, double* fine, + size_t nxc, size_t nyc, size_t nzc, + size_t nxf, size_t nyf, size_t nzf) { + (void)nzc; + for (size_t k = 1; k < nzf - 1; k++) { + size_t K = k >> 1; + size_t nk = (k & 1) ? 2 : 1; + double wk = (k & 1) ? 0.5 : 1.0; + for (size_t j = 1; j < nyf - 1; j++) { + size_t J = j >> 1; + size_t nj = (j & 1) ? 2 : 1; + double wj = (j & 1) ? 0.5 : 1.0; + for (size_t i = 1; i < nxf - 1; i++) { + size_t I = i >> 1; + size_t ni = (i & 1) ? 2 : 1; + double wi = (i & 1) ? 0.5 : 1.0; + double e = 0.0; + for (size_t ck = 0; ck < nk; ck++) { + for (size_t cj = 0; cj < nj; cj++) { + for (size_t ci = 0; ci < ni; ci++) { + e += wi * wj * wk + * coarse[IDX_3D(I + ci, J + cj, K + ck, + nxc, nyc)]; + } + } + } + fine[IDX_3D(i, j, k, nxf, nyf)] += e; + } + } + } +} + +/* ============================================================================ + * INTERIOR MEAN (Neumann nullspace handling) + * ============================================================================ */ + +double mg_interior_mean(const double* f, size_t nx, size_t ny, size_t nz) { + size_t stride_z = (nz > 1) ? nx * ny : 0; + size_t k_start = (nz > 1) ? 1 : 0; + size_t k_end = (nz > 1) ? nz - 1 : 1; + + double sum = 0.0; + size_t count = 0; + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + sum += f[k * stride_z + IDX_2D(i, j, nx)]; + count++; + } + } + } + return (count > 0) ? sum / (double)count : 0.0; +} + +void mg_subtract_interior_mean(double* f, size_t nx, size_t ny, size_t nz) { + double mean = mg_interior_mean(f, nx, ny, nz); + + size_t stride_z = (nz > 1) ? nx * ny : 0; + size_t k_start = (nz > 1) ? 1 : 0; + size_t k_end = (nz > 1) ? nz - 1 : 1; + + for (size_t k = k_start; k < k_end; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + f[k * stride_z + IDX_2D(i, j, nx)] -= mean; + } + } + } +} diff --git a/lib/src/solvers/linear/linear_solver.c b/lib/src/solvers/linear/linear_solver.c index 149817c2..8f8d1776 100644 --- a/lib/src/solvers/linear/linear_solver.c +++ b/lib/src/solvers/linear/linear_solver.c @@ -44,6 +44,13 @@ poisson_solver_params_t poisson_solver_params_default(void) { params.verbose = false; params.preconditioner = POISSON_PRECOND_NONE; params.restart = 0; /* 0 = auto (GMRES_DEFAULT_RESTART); ignored by non-GMRES methods */ + params.mg_cycle = MG_CYCLE_V; + params.mg_smoother = MG_SMOOTHER_REDBLACK_GS; + params.mg_bc = MG_BC_NEUMANN; + params.mg_pre_smooth = 0; /* 0 = default (2) */ + params.mg_post_smooth = 0; /* 0 = default (2) */ + params.mg_coarse_max_iter = 0; /* 0 = default (50) */ + params.mg_max_levels = 0; /* 0 = auto */ return params; } @@ -152,7 +159,9 @@ poisson_solver_t* poisson_solver_create( poisson_solver_method_t method, poisson_solver_backend_t backend) { - /* Auto-select backend if requested */ + /* Auto-select backend if requested (keep the original request: methods + * without a SIMD backend resolve AUTO differently) */ + poisson_solver_backend_t requested = backend; if (backend == POISSON_BACKEND_AUTO) { backend = select_best_backend(); } @@ -253,7 +262,13 @@ poisson_solver_t* poisson_solver_create( } case POISSON_METHOD_MULTIGRID: - /* Not yet implemented */ + /* Only the scalar backend exists. AUTO means "best available", + * which for multigrid IS scalar; explicit SIMD/OMP/GPU requests + * return NULL (no silent fallbacks). */ + if (requested == POISSON_BACKEND_AUTO || + backend == POISSON_BACKEND_SCALAR) { + return create_multigrid_scalar_solver(); + } return NULL; default: @@ -562,6 +577,7 @@ static poisson_solver_t* g_cached_redblack_scalar = NULL; static poisson_solver_t* g_cached_cg_scalar = NULL; static poisson_solver_t* g_cached_cg_omp = NULL; static poisson_solver_t* g_cached_cg_simd = NULL; +static poisson_solver_t* g_cached_mg_scalar = NULL; /** * Cleanup cached solvers (called at program exit) @@ -603,6 +619,10 @@ static void cleanup_cached_solvers(void) { poisson_solver_destroy(g_cached_cg_simd); g_cached_cg_simd = NULL; } + if (g_cached_mg_scalar) { + poisson_solver_destroy(g_cached_mg_scalar); + g_cached_mg_scalar = NULL; + } } int poisson_solve_3d( @@ -670,6 +690,12 @@ int poisson_solve_3d( backend = POISSON_BACKEND_SIMD; break; + case POISSON_SOLVER_MG_SCALAR: + solver_ptr = &g_cached_mg_scalar; + method = POISSON_METHOD_MULTIGRID; + backend = POISSON_BACKEND_SCALAR; + break; + default: CFD_LOG_ERROR("poisson", "poisson_solve_3d: Unknown solver type %d", solver_type); return -1; @@ -696,11 +722,17 @@ int poisson_solve_3d( *solver_ptr = NULL; } - /* Create new solver */ + /* Create new solver. A failed init (e.g. multigrid on non-2^k+1 + * dims) must not leave a broken solver in the cache. */ *solver_ptr = poisson_solver_create(method, backend); if (*solver_ptr) { - poisson_solver_init(*solver_ptr, nx, ny, nz, dx, dy, dz, NULL); + cfd_status_t init_status = + poisson_solver_init(*solver_ptr, nx, ny, nz, dx, dy, dz, NULL); + if (init_status != CFD_SUCCESS) { + poisson_solver_destroy(*solver_ptr); + *solver_ptr = NULL; + } } } diff --git a/lib/src/solvers/linear/linear_solver_internal.h b/lib/src/solvers/linear/linear_solver_internal.h index d9eac5c4..fa4b3f9c 100644 --- a/lib/src/solvers/linear/linear_solver_internal.h +++ b/lib/src/solvers/linear/linear_solver_internal.h @@ -70,6 +70,9 @@ poisson_solver_t* create_gmres_simd_solver(void); poisson_solver_t* create_gmres_omp_solver(void); #endif +/* Geometric multigrid solver (scalar backend only; V/W/F cycles) */ +poisson_solver_t* create_multigrid_scalar_solver(void); + /* ============================================================================ * CG ALGORITHM CONSTANTS * ============================================================================ */ diff --git a/lib/src/solvers/linear/multigrid_internal.h b/lib/src/solvers/linear/multigrid_internal.h new file mode 100644 index 00000000..ec769a38 --- /dev/null +++ b/lib/src/solvers/linear/multigrid_internal.h @@ -0,0 +1,125 @@ +/** + * @file multigrid_internal.h + * @brief Internal data structures and grid-transfer operators for geometric multigrid + * + * Not part of the public API. The transfer operators live in + * cpu/multigrid_transfer.c with internal (non-exported) linkage so the + * operator unit test can compile that translation unit directly. + */ + +#ifndef CFD_MULTIGRID_INTERNAL_H +#define CFD_MULTIGRID_INTERNAL_H + +#include "cfd/solvers/poisson_solver.h" + +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +/* Defaults resolved at init when the corresponding param is 0 */ +#define MG_DEFAULT_PRE_SMOOTH 2 +#define MG_DEFAULT_POST_SMOOTH 2 +#define MG_DEFAULT_COARSE_MAX_ITER 50 + +/* Coarsest-grid size floor per BC mode: a 3x3 Neumann grid's single interior + * point has an identically zero mirror-BC operator (all stencil neighbors + * mirror the center), so Neumann coarsening stops at 5. A 3x3 Dirichlet grid + * is solved exactly by one Gauss-Seidel sweep. */ +#define MG_MIN_COARSE_DIM_DIRICHLET 3 +#define MG_MIN_COARSE_DIM_NEUMANN 5 + +/** + * Per-level grid data for the multigrid hierarchy. + * + * Level 0 (finest) borrows the caller's x/rhs arrays, so its x and rhs + * pointers stay NULL; coarser levels own their buffers. + */ +typedef struct { + size_t nx, ny, nz; /* Grid dimensions at this level (nz==1 for 2D) */ + size_t total; /* nx * ny * nz */ + size_t stride_z; /* nx*ny for 3D, 0 for 2D (branch-free stencils) */ + size_t k_start, k_end; /* z loop bounds ([1,nz-1) for 3D, [0,1) for 2D) */ + double dx2, dy2; /* Squared grid spacings at this level */ + double inv_dz2; /* 1/dz^2 at this level (0.0 for 2D) */ + double inv_factor; /* 1/(2/dx2 + 2/dy2 + 2*inv_dz2): smoother diagonal inverse */ + double* x; /* Correction vector (levels >= 1; NULL at level 0) */ + double* x_temp; /* Jacobi smoother buffer (NULL unless Jacobi smoother) */ + double* rhs; /* Restricted residual / restricted b (levels >= 1; NULL at level 0) */ + double* residual; /* Residual buffer (all levels except coarsest) */ +} mg_level_t; + +/** + * Multigrid solver context (stored in poisson_solver_t.context) + */ +typedef struct { + int num_levels; /* >= 1; [0] = finest, [num_levels-1] = coarsest */ + mg_level_t* levels; + + /* Resolved parameters (defaults already applied) */ + mg_cycle_type_t cycle_type; + mg_smoother_type_t smoother_type; + mg_bc_type_t bc_mode; + int nu1, nu2; /* Pre/post smoothing sweeps */ + int coarse_max_iter; /* Smoother sweeps on the coarsest grid */ + + int fmg_pending; /* 1: next iterate() performs the FMG nested iteration */ + int initialized; +} mg_context_t; + +/** n is of the form 2^k+1 (k>=1): n-1 is a power of two */ +static inline bool mg_is_pow2_plus1(size_t n) { + return (n >= 3) && (((n - 1) & (n - 2)) == 0); +} + +/* ============================================================================ + * GRID-TRANSFER OPERATORS (cpu/multigrid_transfer.c) + * + * Internal linkage across the library; intentionally NOT CFD_LIBRARY_EXPORT. + * All write coarse/fine INTERIOR points only. Restriction never reads fine + * boundary values (coarse-interior stencils stay within the fine interior); + * prolongation reads coarse boundary values, which the caller must have set + * consistently with the BC mode (zero for Dirichlet corrections, mirrored + * for Neumann). + * + * fold_neumann: with zero-gradient BCs the boundary vertex is slaved to its + * interior neighbor (mirror plane at the half-cell), so prolongation reads + * mirrored ghosts. The adjoint of that folded prolongation doubles the + * toward-boundary weight at boundary-adjacent coarse points. Pass 1 in + * Neumann mode to keep restriction the exact adjoint (up to the standard + * 1/4 (2D) or 1/8 (3D) factor); pass 0 for Dirichlet. + * ============================================================================ */ + +/** Full-weighting restriction, 2D (4-2-1)/16 stencil; coarse interior only */ +void mg_restrict_2d(const double* fine, double* coarse, + size_t nxf, size_t nyf, size_t nxc, size_t nyc, + int fold_neumann); + +/** Full-weighting restriction, 3D 27-point (8-4-2-1)/64 stencil */ +void mg_restrict_3d(const double* fine, double* coarse, + size_t nxf, size_t nyf, size_t nzf, + size_t nxc, size_t nyc, size_t nzc, + int fold_neumann); + +/** Bilinear prolongation, adds into fine interior points only */ +void mg_prolongate_add_2d(const double* coarse, double* fine, + size_t nxc, size_t nyc, size_t nxf, size_t nyf); + +/** Trilinear prolongation, adds into fine interior points only */ +void mg_prolongate_add_3d(const double* coarse, double* fine, + size_t nxc, size_t nyc, size_t nzc, + size_t nxf, size_t nyf, size_t nzf); + +/** Mean of f over interior points */ +double mg_interior_mean(const double* f, size_t nx, size_t ny, size_t nz); + +/** Subtract the interior mean from interior points (boundary untouched) */ +void mg_subtract_interior_mean(double* f, size_t nx, size_t ny, size_t nz); + +#ifdef __cplusplus +} +#endif + +#endif /* CFD_MULTIGRID_INTERNAL_H */ diff --git a/tests/math/test_multigrid_convergence.c b/tests/math/test_multigrid_convergence.c new file mode 100644 index 00000000..fa313901 --- /dev/null +++ b/tests/math/test_multigrid_convergence.c @@ -0,0 +1,783 @@ +/** + * @file test_multigrid_convergence.c + * @brief Geometric multigrid solver convergence and API tests + * + * Tests cover: + * - Creation / initialization / backend rules (AUTO -> scalar, others NULL) + * - Grid dimension validation (2^k+1 per active dimension) + * - V(2,2) per-cycle convergence factor < 0.15 (Dirichlet, acceptance) + * - Grid-size-independent cycle counts, 9x9 .. 129x129 (acceptance) + * - Neumann mode agrees with RB-SOR after mean subtraction (nullspace) + * - W-cycle, F-cycle (FMG discretization accuracy), Jacobi smoother + * - 3D grids, mixed dimensions, inhomogeneous Dirichlet data + * - mg_max_levels cap (two-grid method) + * - poisson_solve_3d convenience preset incl. invalid-dims caching + */ + +#include "unity.h" +#include "cfd/solvers/poisson_solver.h" +#include "cfd/core/memory.h" +#include "cfd/core/indexing.h" +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +void setUp(void) {} +void tearDown(void) {} + +/* ============================================================================ + * HELPERS + * ============================================================================ */ + +static double* create_field(size_t total) { + double* field = (double*)cfd_calloc(total, sizeof(double)); + TEST_ASSERT_NOT_NULL_MESSAGE(field, "field allocation failed"); + return field; +} + +/** RHS for -analytic- p = sin(pi x)sin(pi y): f = -2 pi^2 p (p = 0 boundary) */ +static void init_dirichlet_rhs_2d(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + for (size_t j = 0; j < ny; j++) { + double y = (double)j * dy; + for (size_t i = 0; i < nx; i++) { + double x = (double)i * dx; + rhs[IDX_2D(i, j, nx)] = + -2.0 * M_PI * M_PI * sin(M_PI * x) * sin(M_PI * y); + } + } +} + +/** Neumann-compatible RHS f = cos(2 pi x)cos(2 pi y), interior mean removed */ +static void init_neumann_rhs_2d(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + for (size_t j = 0; j < ny; j++) { + double y = (double)j * dy; + for (size_t i = 0; i < nx; i++) { + double x = (double)i * dx; + rhs[IDX_2D(i, j, nx)] = + cos(2.0 * M_PI * x) * cos(2.0 * M_PI * y); + } + } + + double mean = 0.0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + mean += rhs[IDX_2D(i, j, nx)]; + } + } + mean /= (double)((nx - 2) * (ny - 2)); + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + rhs[IDX_2D(i, j, nx)] -= mean; + } + } +} + +static void subtract_interior_mean_2d(double* f, size_t nx, size_t ny) { + double mean = 0.0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + mean += f[IDX_2D(i, j, nx)]; + } + } + mean /= (double)((nx - 2) * (ny - 2)); + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + f[IDX_2D(i, j, nx)] -= mean; + } + } +} + +static void subtract_interior_mean_3d(double* f, size_t nx, size_t ny, + size_t nz) { + double mean = 0.0; + for (size_t k = 1; k < nz - 1; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + mean += f[IDX_3D(i, j, k, nx, ny)]; + } + } + } + mean /= (double)((nx - 2) * (ny - 2) * (nz - 2)); + for (size_t k = 1; k < nz - 1; k++) { + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + f[IDX_3D(i, j, k, nx, ny)] -= mean; + } + } + } +} + +/** + * Create + init a multigrid solver, solve, destroy. Asserts creation and + * init succeed; returns solve status via *stats. + */ +static cfd_status_t solve_mg_2d(double* x, const double* rhs, + size_t nx, size_t ny, double dx, double dy, + const poisson_solver_params_t* params, + poisson_solver_stats_t* stats) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL_MESSAGE(solver, "multigrid create failed"); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + solver, nx, ny, 1, dx, dy, 0.0, params)); + + double* x_temp = create_field(nx * ny); + cfd_status_t status = poisson_solver_solve(solver, x, x_temp, rhs, stats); + cfd_free(x_temp); + poisson_solver_destroy(solver); + return status; +} + +/* ============================================================================ + * CREATION / VALIDATION + * ============================================================================ */ + +void test_mg_create_init_metadata(void) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL_STRING("multigrid_scalar", solver->name); + TEST_ASSERT_EQUAL(POISSON_METHOD_MULTIGRID, solver->method); + TEST_ASSERT_EQUAL(POISSON_BACKEND_SCALAR, solver->backend); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + solver, 33, 33, 1, 1.0 / 32.0, 1.0 / 32.0, 0.0, NULL)); + poisson_solver_destroy(solver); + + /* AUTO resolves to the only (hence best) available backend: scalar */ + solver = poisson_solver_create(POISSON_METHOD_MULTIGRID, + POISSON_BACKEND_AUTO); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(POISSON_BACKEND_SCALAR, solver->backend); + poisson_solver_destroy(solver); + + /* Explicit non-scalar backends: no silent fallbacks */ + TEST_ASSERT_NULL(poisson_solver_create(POISSON_METHOD_MULTIGRID, + POISSON_BACKEND_SIMD)); + TEST_ASSERT_NULL(poisson_solver_create(POISSON_METHOD_MULTIGRID, + POISSON_BACKEND_OMP)); + TEST_ASSERT_NULL(poisson_solver_create(POISSON_METHOD_MULTIGRID, + POISSON_BACKEND_GPU)); +} + +void test_mg_invalid_dims_rejected(void) { + static const size_t bad_dims[][3] = { + {10, 10, 1}, {33, 34, 1}, {32, 32, 1}, {17, 17, 8}, + }; + static const size_t good_dims[][3] = { + {17, 33, 1}, {17, 17, 9}, + }; + + for (size_t n = 0; n < sizeof(bad_dims) / sizeof(bad_dims[0]); n++) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(CFD_ERROR_INVALID, poisson_solver_init( + solver, bad_dims[n][0], bad_dims[n][1], bad_dims[n][2], + 0.1, 0.1, (bad_dims[n][2] > 1) ? 0.1 : 0.0, NULL)); + poisson_solver_destroy(solver); + } + + for (size_t n = 0; n < sizeof(good_dims) / sizeof(good_dims[0]); n++) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + solver, good_dims[n][0], good_dims[n][1], good_dims[n][2], + 0.1, 0.1, (good_dims[n][2] > 1) ? 0.1 : 0.0, NULL)); + poisson_solver_destroy(solver); + } +} + +/* ============================================================================ + * CONVERGENCE FACTOR AND GRID INDEPENDENCE (acceptance criteria) + * ============================================================================ */ + +void test_mg_vcycle_convergence_factor_dirichlet(void) { + const size_t NX = 65, NY = 65; + const double DX = 1.0 / (double)(NX - 1); + const double DY = 1.0 / (double)(NY - 1); + const int NUM_CYCLES = 8; + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + solver, NX, NY, 1, DX, DY, 0.0, ¶ms)); + + double* x = create_field(NX * NY); + double* rhs = create_field(NX * NY); + init_dirichlet_rhs_2d(rhs, NX, NY, DX, DY); + + double res[8]; + for (int n = 0; n < NUM_CYCLES; n++) { + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_iterate( + solver, x, NULL, rhs, &res[n])); + } + + /* Every per-cycle residual reduction after the first cycle must beat the + * V(2,2) Red-Black GS acceptance factor of 0.15 (grid-size-independent) */ + for (int n = 1; n < NUM_CYCLES - 1; n++) { + TEST_ASSERT_TRUE_MESSAGE(res[n] > 0.0, "residual vanished early"); + double factor = res[n + 1] / res[n]; + TEST_ASSERT_TRUE_MESSAGE(factor < 0.15, + "V(2,2) convergence factor >= 0.15"); + } + + cfd_free(x); + cfd_free(rhs); + poisson_solver_destroy(solver); +} + +void test_mg_grid_independence_dirichlet(void) { + static const size_t sizes[] = {9, 17, 33, 65, 129}; + const size_t num_sizes = sizeof(sizes) / sizeof(sizes[0]); + int min_iters = 1000, max_iters = 0; + + for (size_t n = 0; n < num_sizes; n++) { + size_t nx = sizes[n]; + double dx = 1.0 / (double)(nx - 1); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + params.tolerance = 1e-6; + + double* x = create_field(nx * nx); + double* rhs = create_field(nx * nx); + init_dirichlet_rhs_2d(rhs, nx, nx, dx, dx); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, nx, nx, dx, dx, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 15, + "V-cycle count grew beyond 15"); + + if (stats.iterations < min_iters) min_iters = stats.iterations; + if (stats.iterations > max_iters) max_iters = stats.iterations; + + cfd_free(x); + cfd_free(rhs); + } + + /* The key multigrid property: cycle count does not grow with grid size */ + TEST_ASSERT_TRUE_MESSAGE(max_iters - min_iters <= 3, + "cycle count varies with grid size"); +} + +void test_mg_neumann_grid_independence(void) { + static const size_t sizes[] = {17, 33, 65}; + + for (size_t n = 0; n < sizeof(sizes) / sizeof(sizes[0]); n++) { + size_t nx = sizes[n]; + double dx = 1.0 / (double)(nx - 1); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = 1e-6; /* default Neumann mode */ + + double* x = create_field(nx * nx); + double* rhs = create_field(nx * nx); + init_neumann_rhs_2d(rhs, nx, nx, dx, dx); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, nx, nx, dx, dx, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 15, + "Neumann V-cycle count grew beyond 15"); + + cfd_free(x); + cfd_free(rhs); + } +} + +/* ============================================================================ + * NEUMANN MODE VS CG (nullspace-aware comparison) + * ============================================================================ */ + +void test_mg_neumann_matches_rbsor(void) { + /* Reference is RB-SOR, not CG: RB-SOR applies the zero-gradient BC after + * every sweep and so converges to the same mirror-BC discrete system as + * MG. CG freezes boundary values during its Krylov iterations (applying + * the BC only at start/end), which is a different discrete problem. */ + const size_t NX = 33, NY = 33; + const double DX = 1.0 / (double)(NX - 1); + const double DY = 1.0 / (double)(NY - 1); + + double* rhs = create_field(NX * NY); + init_neumann_rhs_2d(rhs, NX, NY, DX, DY); + + poisson_solver_params_t sor_params = poisson_solver_params_default(); + sor_params.tolerance = 1e-10; + sor_params.absolute_tolerance = 1e-12; + sor_params.max_iterations = 20000; + + poisson_solver_t* sor = poisson_solver_create( + POISSON_METHOD_REDBLACK_SOR, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(sor); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + sor, NX, NY, 1, DX, DY, 0.0, &sor_params)); + + double* x_sor = create_field(NX * NY); + double* x_temp = create_field(NX * NY); + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_solve( + sor, x_sor, x_temp, rhs, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + poisson_solver_destroy(sor); + + /* Multigrid, default Neumann mode */ + poisson_solver_params_t mg_params = poisson_solver_params_default(); + mg_params.tolerance = 1e-10; + mg_params.absolute_tolerance = 1e-12; + + double* x_mg = create_field(NX * NY); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x_mg, rhs, NX, NY, DX, DY, &mg_params, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + /* Both solve the same singular system: solutions differ by a constant, + * so compare after removing the interior mean */ + subtract_interior_mean_2d(x_sor, NX, NY); + subtract_interior_mean_2d(x_mg, NX, NY); + + double max_diff = 0.0; + for (size_t j = 1; j < NY - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + double diff = fabs(x_sor[IDX_2D(i, j, NX)] + - x_mg[IDX_2D(i, j, NX)]); + if (diff > max_diff) max_diff = diff; + } + } + TEST_ASSERT_TRUE_MESSAGE(max_diff < 1e-6, + "MG and RB-SOR disagree beyond 1e-6 after mean removal"); + + cfd_free(rhs); + cfd_free(x_sor); + cfd_free(x_mg); + cfd_free(x_temp); +} + +/* ============================================================================ + * CYCLE VARIANTS AND SMOOTHERS + * ============================================================================ */ + +void test_mg_wcycle_converges(void) { + const size_t NX = 65; + const double DX = 1.0 / (double)(NX - 1); + + double* rhs = create_field(NX * NX); + init_dirichlet_rhs_2d(rhs, NX, NX, DX, DX); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + params.tolerance = 1e-6; + + /* V-cycle baseline */ + double* x = create_field(NX * NX); + poisson_solver_stats_t v_stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, NX, NX, DX, DX, ¶ms, &v_stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, v_stats.status); + cfd_free(x); + + /* W-cycle: at least as strong per cycle */ + params.mg_cycle = MG_CYCLE_W; + x = create_field(NX * NX); + poisson_solver_stats_t w_stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, NX, NX, DX, DX, ¶ms, &w_stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, w_stats.status); + TEST_ASSERT_TRUE_MESSAGE(w_stats.iterations <= v_stats.iterations, + "W-cycle needed more cycles than V-cycle"); + cfd_free(x); + + cfd_free(rhs); +} + +void test_mg_fmg_discretization_accuracy(void) { + const size_t NX = 65; + const double DX = 1.0 / (double)(NX - 1); + + double* rhs = create_field(NX * NX); + init_dirichlet_rhs_2d(rhs, NX, NX, DX, DX); + + /* Reference: fully converged V-cycle solve -> discretization error */ + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + params.tolerance = 1e-10; + params.absolute_tolerance = 1e-12; + + double* x_ref = create_field(NX * NX); + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x_ref, rhs, NX, NX, DX, DX, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + double err_ref = 0.0; + for (size_t j = 1; j < NX - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + double exact = sin(M_PI * (double)i * DX) + * sin(M_PI * (double)j * DX); + double err = fabs(x_ref[IDX_2D(i, j, NX)] - exact); + if (err > err_ref) err_ref = err; + } + } + cfd_free(x_ref); + + /* One FMG pass (a single iterate on a fresh F-cycle solver) must land + * within a small factor of the discretization error */ + params.mg_cycle = MG_CYCLE_F; + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + solver, NX, NX, 1, DX, DX, 0.0, ¶ms)); + + double* x_fmg = create_field(NX * NX); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_iterate( + solver, x_fmg, NULL, rhs, NULL)); + + double err_fmg = 0.0; + for (size_t j = 1; j < NX - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + double exact = sin(M_PI * (double)i * DX) + * sin(M_PI * (double)j * DX); + double err = fabs(x_fmg[IDX_2D(i, j, NX)] - exact); + if (err > err_fmg) err_fmg = err; + } + } + TEST_ASSERT_TRUE_MESSAGE(err_fmg <= 2.0 * err_ref, + "one FMG pass missed discretization accuracy"); + cfd_free(x_fmg); + poisson_solver_destroy(solver); + + /* A full F-cycle solve converges in very few cycles */ + params.tolerance = 1e-6; + params.absolute_tolerance = 1e-10; + double* x = create_field(NX * NX); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, NX, NX, DX, DX, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 5, + "F-cycle solve took more than 5 cycles"); + cfd_free(x); + + cfd_free(rhs); +} + +void test_mg_jacobi_smoother_converges(void) { + const size_t NX = 33; + const double DX = 1.0 / (double)(NX - 1); + + double* x = create_field(NX * NX); + double* rhs = create_field(NX * NX); + init_dirichlet_rhs_2d(rhs, NX, NX, DX, DX); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + params.mg_smoother = MG_SMOOTHER_JACOBI; + params.tolerance = 1e-6; + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, NX, NX, DX, DX, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + /* Weighted Jacobi smooths less per sweep than RB-GS: looser bound */ + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 25, + "Jacobi-smoothed V-cycle count grew beyond 25"); + + cfd_free(x); + cfd_free(rhs); +} + +/* ============================================================================ + * 3D AND MIXED DIMENSIONS + * ============================================================================ */ + +static void run_mg_3d_case(size_t n, mg_bc_type_t bc, int compare_rbsor) { + double d = 1.0 / (double)(n - 1); + size_t total = n * n * n; + + double* rhs = create_field(total); + for (size_t k = 0; k < n; k++) { + double z = (double)k * d; + for (size_t j = 0; j < n; j++) { + double y = (double)j * d; + for (size_t i = 0; i < n; i++) { + double x = (double)i * d; + if (bc == MG_BC_DIRICHLET) { + rhs[IDX_3D(i, j, k, n, n)] = -3.0 * M_PI * M_PI + * sin(M_PI * x) * sin(M_PI * y) * sin(M_PI * z); + } else { + rhs[IDX_3D(i, j, k, n, n)] = cos(2.0 * M_PI * x) + * cos(2.0 * M_PI * y) * cos(2.0 * M_PI * z); + } + } + } + } + if (bc == MG_BC_NEUMANN) { + subtract_interior_mean_3d(rhs, n, n, n); + } + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = bc; + params.tolerance = 1e-6; + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + solver, n, n, n, d, d, d, ¶ms)); + + double* x = create_field(total); + double* x_temp = create_field(total); + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_solve( + solver, x, x_temp, rhs, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 15, + "3D V-cycle count grew beyond 15"); + poisson_solver_destroy(solver); + + if (compare_rbsor) { + /* RB-SOR solves the same mirror-BC system as MG (see 2D test) */ + poisson_solver_params_t sor_params = poisson_solver_params_default(); + sor_params.tolerance = 1e-10; + sor_params.absolute_tolerance = 1e-12; + sor_params.max_iterations = 20000; + + poisson_solver_t* sor = poisson_solver_create( + POISSON_METHOD_REDBLACK_SOR, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(sor); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_init( + sor, n, n, n, d, d, d, &sor_params)); + + double* x_sor = create_field(total); + TEST_ASSERT_EQUAL(CFD_SUCCESS, poisson_solver_solve( + sor, x_sor, x_temp, rhs, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + poisson_solver_destroy(sor); + + subtract_interior_mean_3d(x, n, n, n); + subtract_interior_mean_3d(x_sor, n, n, n); + + double max_diff = 0.0; + for (size_t k = 1; k < n - 1; k++) { + for (size_t j = 1; j < n - 1; j++) { + for (size_t i = 1; i < n - 1; i++) { + double diff = fabs(x[IDX_3D(i, j, k, n, n)] + - x_sor[IDX_3D(i, j, k, n, n)]); + if (diff > max_diff) max_diff = diff; + } + } + } + TEST_ASSERT_TRUE_MESSAGE(max_diff < 1e-5, + "3D MG and RB-SOR disagree after mean removal"); + cfd_free(x_sor); + } + + cfd_free(x); + cfd_free(x_temp); + cfd_free(rhs); +} + +void test_mg_3d_convergence(void) { + run_mg_3d_case(9, MG_BC_DIRICHLET, 0); + run_mg_3d_case(17, MG_BC_DIRICHLET, 0); + run_mg_3d_case(9, MG_BC_NEUMANN, 0); + run_mg_3d_case(17, MG_BC_NEUMANN, 1); +} + +void test_mg_mixed_dims(void) { + static const size_t dims[][2] = {{17, 33}, {33, 17}}; + + for (size_t n = 0; n < sizeof(dims) / sizeof(dims[0]); n++) { + size_t nx = dims[n][0]; + size_t ny = dims[n][1]; + double dx = 1.0 / (double)(nx - 1); + double dy = 1.0 / (double)(ny - 1); + + for (int bc = 0; bc < 2; bc++) { + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = (bc == 0) ? MG_BC_NEUMANN : MG_BC_DIRICHLET; + params.tolerance = 1e-6; + + double* x = create_field(nx * ny); + double* rhs = create_field(nx * ny); + if (params.mg_bc == MG_BC_DIRICHLET) { + init_dirichlet_rhs_2d(rhs, nx, ny, dx, dy); + } else { + init_neumann_rhs_2d(rhs, nx, ny, dx, dy); + } + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, nx, ny, dx, dy, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 15, + "mixed-dims cycle count grew beyond 15"); + + cfd_free(x); + cfd_free(rhs); + } + } +} + +/* ============================================================================ + * INHOMOGENEOUS DIRICHLET DATA + * ============================================================================ */ + +void test_mg_dirichlet_inhomogeneous(void) { + /* p = x^2 + y^2 has Laplacian(p) = 4 represented EXACTLY by the 5-point + * stencil, so the converged discrete solution equals the analytic one. + * Verifies caller-supplied boundary values are held fixed and enter the + * interior stencil correctly. */ + const size_t NX = 17; + const double DX = 1.0 / (double)(NX - 1); + + double* x = create_field(NX * NX); + double* rhs = create_field(NX * NX); + for (size_t j = 0; j < NX; j++) { + double yv = (double)j * DX; + for (size_t i = 0; i < NX; i++) { + double xv = (double)i * DX; + rhs[IDX_2D(i, j, NX)] = 4.0; + if (i == 0 || j == 0 || i == NX - 1 || j == NX - 1) { + x[IDX_2D(i, j, NX)] = xv * xv + yv * yv; + } + } + } + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + params.tolerance = 1e-12; + params.absolute_tolerance = 1e-12; + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, NX, NX, DX, DX, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + + for (size_t j = 1; j < NX - 1; j++) { + double yv = (double)j * DX; + for (size_t i = 1; i < NX - 1; i++) { + double xv = (double)i * DX; + TEST_ASSERT_DOUBLE_WITHIN(1e-8, xv * xv + yv * yv, + x[IDX_2D(i, j, NX)]); + } + } + + cfd_free(x); + cfd_free(rhs); +} + +/* ============================================================================ + * LEVEL CAP AND CONVENIENCE API + * ============================================================================ */ + +void test_mg_max_levels_cap(void) { + const size_t NX = 65; + const double DX = 1.0 / (double)(NX - 1); + + double* x = create_field(NX * NX); + double* rhs = create_field(NX * NX); + init_dirichlet_rhs_2d(rhs, NX, NX, DX, DX); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.mg_bc = MG_BC_DIRICHLET; + params.mg_max_levels = 2; + /* Two-grid: the 33x33 "coarsest" level needs a real solve, not 50 sweeps */ + params.mg_coarse_max_iter = 500; + params.tolerance = 1e-6; + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + TEST_ASSERT_EQUAL(CFD_SUCCESS, solve_mg_2d( + x, rhs, NX, NX, DX, DX, ¶ms, &stats)); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats.status); + TEST_ASSERT_TRUE_MESSAGE(stats.iterations <= 60, + "two-grid method failed to converge in 60 cycles"); + + cfd_free(x); + cfd_free(rhs); +} + +void test_mg_convenience_api(void) { + const size_t NX = 33; + const double DX = 1.0 / (double)(NX - 1); + + double* p = create_field(NX * NX); + double* p_temp = create_field(NX * NX); + double* rhs = create_field(NX * NX); + init_neumann_rhs_2d(rhs, NX, NX, DX, DX); + + int iters = poisson_solve_3d(p, p_temp, rhs, NX, NX, 1, DX, DX, 0.0, + POISSON_SOLVER_MG_SCALAR); + TEST_ASSERT_TRUE_MESSAGE(iters > 0, "convenience MG solve failed"); + + /* Second call with same dims exercises the cached-solver path */ + for (size_t n = 0; n < NX * NX; n++) { + p[n] = 0.0; + } + iters = poisson_solve_3d(p, p_temp, rhs, NX, NX, 1, DX, DX, 0.0, + POISSON_SOLVER_MG_SCALAR); + TEST_ASSERT_TRUE_MESSAGE(iters > 0, "cached convenience MG solve failed"); + + /* Invalid dims: init fails, cache slot must not keep a broken solver */ + double* p32 = create_field(32 * 32); + double* rhs32 = create_field(32 * 32); + iters = poisson_solve_3d(p32, NULL, rhs32, 32, 32, 1, + 1.0 / 31.0, 1.0 / 31.0, 0.0, + POISSON_SOLVER_MG_SCALAR); + TEST_ASSERT_EQUAL(-1, iters); + /* And a valid solve afterwards still works (cache not poisoned) */ + for (size_t n = 0; n < NX * NX; n++) { + p[n] = 0.0; + } + iters = poisson_solve_3d(p, p_temp, rhs, NX, NX, 1, DX, DX, 0.0, + POISSON_SOLVER_MG_SCALAR); + TEST_ASSERT_TRUE_MESSAGE(iters > 0, "MG solve after failed init broke"); + + cfd_free(p); + cfd_free(p_temp); + cfd_free(rhs); + cfd_free(p32); + cfd_free(rhs32); +} + +/* ============================================================================ + * MAIN + * ============================================================================ */ + +int main(void) { + UNITY_BEGIN(); + + RUN_TEST(test_mg_create_init_metadata); + RUN_TEST(test_mg_invalid_dims_rejected); + + RUN_TEST(test_mg_vcycle_convergence_factor_dirichlet); + RUN_TEST(test_mg_grid_independence_dirichlet); + RUN_TEST(test_mg_neumann_grid_independence); + RUN_TEST(test_mg_neumann_matches_rbsor); + + RUN_TEST(test_mg_wcycle_converges); + RUN_TEST(test_mg_fmg_discretization_accuracy); + RUN_TEST(test_mg_jacobi_smoother_converges); + + RUN_TEST(test_mg_3d_convergence); + RUN_TEST(test_mg_mixed_dims); + RUN_TEST(test_mg_dirichlet_inhomogeneous); + + RUN_TEST(test_mg_max_levels_cap); + RUN_TEST(test_mg_convenience_api); + + return UNITY_END(); +} diff --git a/tests/math/test_multigrid_operators.c b/tests/math/test_multigrid_operators.c new file mode 100644 index 00000000..5114d2d3 --- /dev/null +++ b/tests/math/test_multigrid_operators.c @@ -0,0 +1,681 @@ +/** + * @file test_multigrid_operators.c + * @brief Unit tests for multigrid grid-transfer operators + * + * Compiles lib/src/solvers/linear/cpu/multigrid_transfer.c directly (the + * operators have internal linkage and are not exported from the library). + * Deliberately uses no poisson_solver_* API and no cfd runtime — pure math + * verification of the operators before they are used in the solver. + * + * Tests cover: + * - Restriction exactness for constant and linear fields (2D/3D) + * - Restriction writes coarse interior only (boundary contract) + * - Prolongation exactness for constant and linear fields (2D/3D) + * - Prolongation adds into fine interior only (boundary contract) + * - Adjointness: = (1/4) in 2D, (1/8) in 3D + * - Conservation: sum(Rf) = (1/4)sum(f) in 2D, (1/8) in 3D + * - Interior mean subtraction + * - 2^k+1 dimension predicate + */ + +#include "unity.h" +#include "../../lib/src/solvers/linear/multigrid_internal.h" +#include "cfd/core/indexing.h" + +#include +#include + +void setUp(void) {} +void tearDown(void) {} + +#define SENTINEL 999.0 + +/* ============================================================================ + * HELPERS + * ============================================================================ */ + +/** Deterministic LCG in [-1, 1] (no rand() — reproducible across platforms) */ +static unsigned int g_lcg_state = 12345u; + +static double lcg_next(void) { + g_lcg_state = g_lcg_state * 1664525u + 1013904223u; + return (double)(g_lcg_state >> 8) / (double)(1u << 24) * 2.0 - 1.0; +} + +static double* alloc_zeroed(size_t n) { + double* p = (double*)calloc(n, sizeof(double)); + TEST_ASSERT_NOT_NULL_MESSAGE(p, "allocation failed"); + return p; +} + +/** Fill a 2D field with f(x,y) = a + b*x + c*y on the unit square */ +static void fill_linear_2d(double* f, size_t nx, size_t ny, + double a, double b, double c) { + double dx = 1.0 / (double)(nx - 1); + double dy = 1.0 / (double)(ny - 1); + for (size_t j = 0; j < ny; j++) { + for (size_t i = 0; i < nx; i++) { + f[IDX_2D(i, j, nx)] = a + b * (double)i * dx + c * (double)j * dy; + } + } +} + +/** Fill a 3D field with f(x,y,z) = a + b*x + c*y + d*z on the unit cube */ +static void fill_linear_3d(double* f, size_t nx, size_t ny, size_t nz, + double a, double b, double c, double d) { + double dx = 1.0 / (double)(nx - 1); + double dy = 1.0 / (double)(ny - 1); + double dz = 1.0 / (double)(nz - 1); + for (size_t k = 0; k < nz; k++) { + for (size_t j = 0; j < ny; j++) { + for (size_t i = 0; i < nx; i++) { + f[IDX_3D(i, j, k, nx, ny)] = a + b * (double)i * dx + + c * (double)j * dy + + d * (double)k * dz; + } + } + } +} + +static int is_boundary_2d(size_t i, size_t j, size_t nx, size_t ny) { + return i == 0 || j == 0 || i == nx - 1 || j == ny - 1; +} + +static int is_boundary_3d(size_t i, size_t j, size_t k, + size_t nx, size_t ny, size_t nz) { + return is_boundary_2d(i, j, nx, ny) || k == 0 || k == nz - 1; +} + +/* ============================================================================ + * RESTRICTION TESTS (2D) + * ============================================================================ */ + +void test_restrict_2d_constant_exact(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* fine = alloc_zeroed(NXF * NYF); + double* coarse = alloc_zeroed(NXC * NYC); + + fill_linear_2d(fine, NXF, NYF, 3.7, 0.0, 0.0); + mg_restrict_2d(fine, coarse, NXF, NYF, NXC, NYC, 0); + + for (size_t J = 1; J < NYC - 1; J++) { + for (size_t I = 1; I < NXC - 1; I++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-14, 3.7, coarse[IDX_2D(I, J, NXC)]); + } + } + + free(fine); + free(coarse); +} + +void test_restrict_2d_linear_exact(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* fine = alloc_zeroed(NXF * NYF); + double* coarse = alloc_zeroed(NXC * NYC); + + fill_linear_2d(fine, NXF, NYF, 1.0, 2.0, 3.0); + mg_restrict_2d(fine, coarse, NXF, NYF, NXC, NYC, 0); + + /* Coarse node (I,J) is coincident with fine node (2I,2J); full weighting + * of a linear field returns the center value exactly (symmetric stencil) */ + for (size_t J = 1; J < NYC - 1; J++) { + for (size_t I = 1; I < NXC - 1; I++) { + double expected = fine[IDX_2D(2 * I, 2 * J, NXF)]; + TEST_ASSERT_DOUBLE_WITHIN(1e-13, expected, + coarse[IDX_2D(I, J, NXC)]); + } + } + + free(fine); + free(coarse); +} + +void test_restrict_2d_boundary_contract(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* fine = alloc_zeroed(NXF * NYF); + double* coarse = alloc_zeroed(NXC * NYC); + + for (size_t i = 0; i < NXF * NYF; i++) { + fine[i] = lcg_next(); + } + for (size_t J = 0; J < NYC; J++) { + for (size_t I = 0; I < NXC; I++) { + if (is_boundary_2d(I, J, NXC, NYC)) { + coarse[IDX_2D(I, J, NXC)] = SENTINEL; + } + } + } + + mg_restrict_2d(fine, coarse, NXF, NYF, NXC, NYC, 0); + + for (size_t J = 0; J < NYC; J++) { + for (size_t I = 0; I < NXC; I++) { + if (is_boundary_2d(I, J, NXC, NYC)) { + TEST_ASSERT_EQUAL_DOUBLE(SENTINEL, coarse[IDX_2D(I, J, NXC)]); + } + } + } + + free(fine); + free(coarse); +} + +/* ============================================================================ + * RESTRICTION TESTS (3D) + * ============================================================================ */ + +void test_restrict_3d_constant_and_linear_exact(void) { + const size_t NF = 9, NC = 5; + double* fine = alloc_zeroed(NF * NF * NF); + double* coarse = alloc_zeroed(NC * NC * NC); + + /* Constant */ + fill_linear_3d(fine, NF, NF, NF, -2.5, 0.0, 0.0, 0.0); + mg_restrict_3d(fine, coarse, NF, NF, NF, NC, NC, NC, 0); + for (size_t K = 1; K < NC - 1; K++) { + for (size_t J = 1; J < NC - 1; J++) { + for (size_t I = 1; I < NC - 1; I++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-14, -2.5, + coarse[IDX_3D(I, J, K, NC, NC)]); + } + } + } + + /* Linear */ + fill_linear_3d(fine, NF, NF, NF, 1.0, 2.0, 3.0, 4.0); + mg_restrict_3d(fine, coarse, NF, NF, NF, NC, NC, NC, 0); + for (size_t K = 1; K < NC - 1; K++) { + for (size_t J = 1; J < NC - 1; J++) { + for (size_t I = 1; I < NC - 1; I++) { + double expected = fine[IDX_3D(2 * I, 2 * J, 2 * K, NF, NF)]; + TEST_ASSERT_DOUBLE_WITHIN(1e-13, expected, + coarse[IDX_3D(I, J, K, NC, NC)]); + } + } + } + + free(fine); + free(coarse); +} + +/* ============================================================================ + * PROLONGATION TESTS (2D) + * ============================================================================ */ + +void test_prolongate_2d_constant_exact(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* fine = alloc_zeroed(NXF * NYF); + double* coarse = alloc_zeroed(NXC * NYC); + + fill_linear_2d(coarse, NXC, NYC, 4.2, 0.0, 0.0); /* incl. boundary */ + mg_prolongate_add_2d(coarse, fine, NXC, NYC, NXF, NYF); + + for (size_t j = 0; j < NYF; j++) { + for (size_t i = 0; i < NXF; i++) { + if (is_boundary_2d(i, j, NXF, NYF)) { + /* Additive prolongation must not touch the fine boundary */ + TEST_ASSERT_EQUAL_DOUBLE(0.0, fine[IDX_2D(i, j, NXF)]); + } else { + TEST_ASSERT_DOUBLE_WITHIN(1e-14, 4.2, fine[IDX_2D(i, j, NXF)]); + } + } + } + + free(fine); + free(coarse); +} + +void test_prolongate_2d_linear_exact(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* fine = alloc_zeroed(NXF * NYF); + double* coarse = alloc_zeroed(NXC * NYC); + double* expected = alloc_zeroed(NXF * NYF); + + fill_linear_2d(coarse, NXC, NYC, 1.0, 2.0, 3.0); + fill_linear_2d(expected, NXF, NYF, 1.0, 2.0, 3.0); + mg_prolongate_add_2d(coarse, fine, NXC, NYC, NXF, NYF); + + /* Bilinear interpolation is exact for linear functions */ + for (size_t j = 1; j < NYF - 1; j++) { + for (size_t i = 1; i < NXF - 1; i++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-13, expected[IDX_2D(i, j, NXF)], + fine[IDX_2D(i, j, NXF)]); + } + } + + free(fine); + free(coarse); + free(expected); +} + +/* ============================================================================ + * PROLONGATION TESTS (3D) + * ============================================================================ */ + +void test_prolongate_3d_linear_exact(void) { + const size_t NF = 9, NC = 5; + double* fine = alloc_zeroed(NF * NF * NF); + double* coarse = alloc_zeroed(NC * NC * NC); + double* expected = alloc_zeroed(NF * NF * NF); + + fill_linear_3d(coarse, NC, NC, NC, 1.0, 2.0, 3.0, 4.0); + fill_linear_3d(expected, NF, NF, NF, 1.0, 2.0, 3.0, 4.0); + mg_prolongate_add_3d(coarse, fine, NC, NC, NC, NF, NF, NF); + + for (size_t k = 0; k < NF; k++) { + for (size_t j = 0; j < NF; j++) { + for (size_t i = 0; i < NF; i++) { + if (is_boundary_3d(i, j, k, NF, NF, NF)) { + TEST_ASSERT_EQUAL_DOUBLE(0.0, + fine[IDX_3D(i, j, k, NF, NF)]); + } else { + TEST_ASSERT_DOUBLE_WITHIN(1e-13, + expected[IDX_3D(i, j, k, NF, NF)], + fine[IDX_3D(i, j, k, NF, NF)]); + } + } + } + } + + free(fine); + free(coarse); + free(expected); +} + +/* ============================================================================ + * ADJOINTNESS: = c * (c = 1/4 in 2D, 1/8 in 3D) + * + * Essential for the coarse-grid operator to inherit positive-definiteness. + * f, g have zero boundary values so the interior-truncated operators satisfy + * the identity exactly. + * ============================================================================ */ + +void test_adjointness_2d(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* f = alloc_zeroed(NXF * NYF); + double* g = alloc_zeroed(NXC * NYC); + double* Rf = alloc_zeroed(NXC * NYC); + double* Pg = alloc_zeroed(NXF * NYF); + + for (size_t j = 1; j < NYF - 1; j++) { + for (size_t i = 1; i < NXF - 1; i++) { + f[IDX_2D(i, j, NXF)] = lcg_next(); + } + } + for (size_t J = 1; J < NYC - 1; J++) { + for (size_t I = 1; I < NXC - 1; I++) { + g[IDX_2D(I, J, NXC)] = lcg_next(); + } + } + + mg_restrict_2d(f, Rf, NXF, NYF, NXC, NYC, 0); + mg_prolongate_add_2d(g, Pg, NXC, NYC, NXF, NYF); + + double lhs = 0.0, rhs = 0.0; + for (size_t n = 0; n < NXC * NYC; n++) { + lhs += Rf[n] * g[n]; + } + for (size_t n = 0; n < NXF * NYF; n++) { + rhs += f[n] * Pg[n]; + } + rhs *= 0.25; + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(rhs) + 1e-14, rhs, lhs); + + free(f); + free(g); + free(Rf); + free(Pg); +} + +void test_adjointness_3d(void) { + const size_t NF = 9, NC = 5; + double* f = alloc_zeroed(NF * NF * NF); + double* g = alloc_zeroed(NC * NC * NC); + double* Rf = alloc_zeroed(NC * NC * NC); + double* Pg = alloc_zeroed(NF * NF * NF); + + for (size_t k = 1; k < NF - 1; k++) { + for (size_t j = 1; j < NF - 1; j++) { + for (size_t i = 1; i < NF - 1; i++) { + f[IDX_3D(i, j, k, NF, NF)] = lcg_next(); + } + } + } + for (size_t K = 1; K < NC - 1; K++) { + for (size_t J = 1; J < NC - 1; J++) { + for (size_t I = 1; I < NC - 1; I++) { + g[IDX_3D(I, J, K, NC, NC)] = lcg_next(); + } + } + } + + mg_restrict_3d(f, Rf, NF, NF, NF, NC, NC, NC, 0); + mg_prolongate_add_3d(g, Pg, NC, NC, NC, NF, NF, NF); + + double lhs = 0.0, rhs = 0.0; + for (size_t n = 0; n < NC * NC * NC; n++) { + lhs += Rf[n] * g[n]; + } + for (size_t n = 0; n < NF * NF * NF; n++) { + rhs += f[n] * Pg[n]; + } + rhs *= 0.125; + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(rhs) + 1e-14, rhs, lhs); + + free(f); + free(g); + free(Rf); + free(Pg); +} + +/* ============================================================================ + * CONSERVATION: sum(Rf) = c * sum(f) (c = 1/4 in 2D, 1/8 in 3D) + * + * Holds when f is supported away from the first interior ring, so every + * nonzero fine point has its full restriction-stencil coverage (total + * weight exactly 1/4 per fine point in 2D, 1/8 in 3D). + * ============================================================================ */ + +void test_conservation_2d(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* f = alloc_zeroed(NXF * NYF); + double* Rf = alloc_zeroed(NXC * NYC); + + double sum_f = 0.0; + for (size_t j = 2; j < NYF - 2; j++) { + for (size_t i = 2; i < NXF - 2; i++) { + double v = lcg_next(); + f[IDX_2D(i, j, NXF)] = v; + sum_f += v; + } + } + + mg_restrict_2d(f, Rf, NXF, NYF, NXC, NYC, 0); + + double sum_Rf = 0.0; + for (size_t n = 0; n < NXC * NYC; n++) { + sum_Rf += Rf[n]; + } + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(sum_f) + 1e-14, + 0.25 * sum_f, sum_Rf); + + free(f); + free(Rf); +} + +void test_conservation_3d(void) { + const size_t NF = 9, NC = 5; + double* f = alloc_zeroed(NF * NF * NF); + double* Rf = alloc_zeroed(NC * NC * NC); + + double sum_f = 0.0; + for (size_t k = 2; k < NF - 2; k++) { + for (size_t j = 2; j < NF - 2; j++) { + for (size_t i = 2; i < NF - 2; i++) { + double v = lcg_next(); + f[IDX_3D(i, j, k, NF, NF)] = v; + sum_f += v; + } + } + } + + mg_restrict_3d(f, Rf, NF, NF, NF, NC, NC, NC, 0); + + double sum_Rf = 0.0; + for (size_t n = 0; n < NC * NC * NC; n++) { + sum_Rf += Rf[n]; + } + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(sum_f) + 1e-14, + 0.125 * sum_f, sum_Rf); + + free(f); + free(Rf); +} + +/* ============================================================================ + * NEUMANN-FOLDED RESTRICTION (fold_neumann = 1) + * + * With zero-gradient BCs, prolongation reads mirrored ghost values (boundary + * = adjacent interior; corners = diagonal interior). The folded restriction + * must be the exact adjoint of that folded prolongation, and constants must + * be reproduced by the folded prolongation, which makes the conservation + * identity hold for ARBITRARY interior fields. + * ============================================================================ */ + +/** Mirror ghost fill matching bc_apply_scalar Neumann: x-faces then y-faces + * (corners end up as diagonal interior copies) */ +static void mirror_fill_2d(double* f, size_t nx, size_t ny) { + for (size_t j = 0; j < ny; j++) { + f[IDX_2D(0, j, nx)] = f[IDX_2D(1, j, nx)]; + f[IDX_2D(nx - 1, j, nx)] = f[IDX_2D(nx - 2, j, nx)]; + } + for (size_t i = 0; i < nx; i++) { + f[IDX_2D(i, 0, nx)] = f[IDX_2D(i, 1, nx)]; + f[IDX_2D(i, ny - 1, nx)] = f[IDX_2D(i, ny - 2, nx)]; + } +} + +/** 3D mirror fill: x-faces, y-faces per plane, then z-plane copies */ +static void mirror_fill_3d(double* f, size_t nx, size_t ny, size_t nz) { + for (size_t k = 0; k < nz; k++) { + mirror_fill_2d(f + k * nx * ny, nx, ny); + } + for (size_t n = 0; n < nx * ny; n++) { + f[n] = f[nx * ny + n]; + f[(nz - 1) * nx * ny + n] = f[(nz - 2) * nx * ny + n]; + } +} + +void test_adjointness_2d_folded(void) { + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* f = alloc_zeroed(NXF * NYF); + double* g = alloc_zeroed(NXC * NYC); + double* Rf = alloc_zeroed(NXC * NYC); + double* Pg = alloc_zeroed(NXF * NYF); + + for (size_t j = 1; j < NYF - 1; j++) { + for (size_t i = 1; i < NXF - 1; i++) { + f[IDX_2D(i, j, NXF)] = lcg_next(); + } + } + for (size_t J = 1; J < NYC - 1; J++) { + for (size_t I = 1; I < NXC - 1; I++) { + g[IDX_2D(I, J, NXC)] = lcg_next(); + } + } + + mg_restrict_2d(f, Rf, NXF, NYF, NXC, NYC, 1); + mirror_fill_2d(g, NXC, NYC); + mg_prolongate_add_2d(g, Pg, NXC, NYC, NXF, NYF); + + double lhs = 0.0, rhs = 0.0; + for (size_t J = 1; J < NYC - 1; J++) { + for (size_t I = 1; I < NXC - 1; I++) { + lhs += Rf[IDX_2D(I, J, NXC)] * g[IDX_2D(I, J, NXC)]; + } + } + for (size_t j = 1; j < NYF - 1; j++) { + for (size_t i = 1; i < NXF - 1; i++) { + rhs += f[IDX_2D(i, j, NXF)] * Pg[IDX_2D(i, j, NXF)]; + } + } + rhs *= 0.25; + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(rhs) + 1e-14, rhs, lhs); + + free(f); + free(g); + free(Rf); + free(Pg); +} + +void test_adjointness_3d_folded(void) { + const size_t NF = 9, NC = 5; + double* f = alloc_zeroed(NF * NF * NF); + double* g = alloc_zeroed(NC * NC * NC); + double* Rf = alloc_zeroed(NC * NC * NC); + double* Pg = alloc_zeroed(NF * NF * NF); + + for (size_t k = 1; k < NF - 1; k++) { + for (size_t j = 1; j < NF - 1; j++) { + for (size_t i = 1; i < NF - 1; i++) { + f[IDX_3D(i, j, k, NF, NF)] = lcg_next(); + } + } + } + for (size_t K = 1; K < NC - 1; K++) { + for (size_t J = 1; J < NC - 1; J++) { + for (size_t I = 1; I < NC - 1; I++) { + g[IDX_3D(I, J, K, NC, NC)] = lcg_next(); + } + } + } + + mg_restrict_3d(f, Rf, NF, NF, NF, NC, NC, NC, 1); + mirror_fill_3d(g, NC, NC, NC); + mg_prolongate_add_3d(g, Pg, NC, NC, NC, NF, NF, NF); + + double lhs = 0.0, rhs = 0.0; + for (size_t K = 1; K < NC - 1; K++) { + for (size_t J = 1; J < NC - 1; J++) { + for (size_t I = 1; I < NC - 1; I++) { + lhs += Rf[IDX_3D(I, J, K, NC, NC)] + * g[IDX_3D(I, J, K, NC, NC)]; + } + } + } + for (size_t k = 1; k < NF - 1; k++) { + for (size_t j = 1; j < NF - 1; j++) { + for (size_t i = 1; i < NF - 1; i++) { + rhs += f[IDX_3D(i, j, k, NF, NF)] + * Pg[IDX_3D(i, j, k, NF, NF)]; + } + } + } + rhs *= 0.125; + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(rhs) + 1e-14, rhs, lhs); + + free(f); + free(g); + free(Rf); + free(Pg); +} + +void test_conservation_2d_folded(void) { + /* Folded prolongation reproduces constants exactly, so conservation + * holds for ANY interior field (no support restriction) */ + const size_t NXF = 17, NYF = 17, NXC = 9, NYC = 9; + double* f = alloc_zeroed(NXF * NYF); + double* Rf = alloc_zeroed(NXC * NYC); + + double sum_f = 0.0; + for (size_t j = 1; j < NYF - 1; j++) { + for (size_t i = 1; i < NXF - 1; i++) { + double v = lcg_next(); + f[IDX_2D(i, j, NXF)] = v; + sum_f += v; + } + } + + mg_restrict_2d(f, Rf, NXF, NYF, NXC, NYC, 1); + + double sum_Rf = 0.0; + for (size_t J = 1; J < NYC - 1; J++) { + for (size_t I = 1; I < NXC - 1; I++) { + sum_Rf += Rf[IDX_2D(I, J, NXC)]; + } + } + + TEST_ASSERT_DOUBLE_WITHIN(1e-12 * fabs(sum_f) + 1e-14, + 0.25 * sum_f, sum_Rf); + + free(f); + free(Rf); +} + +/* ============================================================================ + * INTERIOR MEAN + * ============================================================================ */ + +void test_subtract_interior_mean(void) { + const size_t NX = 17, NY = 9; + double* f = alloc_zeroed(NX * NY); + + for (size_t n = 0; n < NX * NY; n++) { + f[n] = 5.0 + lcg_next(); + } + for (size_t j = 0; j < NY; j++) { + for (size_t i = 0; i < NX; i++) { + if (is_boundary_2d(i, j, NX, NY)) { + f[IDX_2D(i, j, NX)] = SENTINEL; + } + } + } + + mg_subtract_interior_mean(f, NX, NY, 1); + + TEST_ASSERT_DOUBLE_WITHIN(1e-14, 0.0, mg_interior_mean(f, NX, NY, 1)); + for (size_t j = 0; j < NY; j++) { + for (size_t i = 0; i < NX; i++) { + if (is_boundary_2d(i, j, NX, NY)) { + TEST_ASSERT_EQUAL_DOUBLE(SENTINEL, f[IDX_2D(i, j, NX)]); + } + } + } + + free(f); +} + +/* ============================================================================ + * DIMENSION PREDICATE + * ============================================================================ */ + +void test_is_pow2_plus1(void) { + static const size_t valid[] = {3, 5, 9, 17, 33, 65, 129, 257}; + static const size_t invalid[] = {0, 1, 2, 4, 6, 7, 8, 10, 16, 32, 34}; + + for (size_t n = 0; n < sizeof(valid) / sizeof(valid[0]); n++) { + TEST_ASSERT_TRUE_MESSAGE(mg_is_pow2_plus1(valid[n]), + "expected 2^k+1 dimension to be accepted"); + } + for (size_t n = 0; n < sizeof(invalid) / sizeof(invalid[0]); n++) { + TEST_ASSERT_FALSE_MESSAGE(mg_is_pow2_plus1(invalid[n]), + "expected non-2^k+1 dimension to be rejected"); + } +} + +/* ============================================================================ + * MAIN + * ============================================================================ */ + +int main(void) { + UNITY_BEGIN(); + + RUN_TEST(test_restrict_2d_constant_exact); + RUN_TEST(test_restrict_2d_linear_exact); + RUN_TEST(test_restrict_2d_boundary_contract); + RUN_TEST(test_restrict_3d_constant_and_linear_exact); + + RUN_TEST(test_prolongate_2d_constant_exact); + RUN_TEST(test_prolongate_2d_linear_exact); + RUN_TEST(test_prolongate_3d_linear_exact); + + RUN_TEST(test_adjointness_2d); + RUN_TEST(test_adjointness_3d); + RUN_TEST(test_conservation_2d); + RUN_TEST(test_conservation_3d); + + RUN_TEST(test_adjointness_2d_folded); + RUN_TEST(test_adjointness_3d_folded); + RUN_TEST(test_conservation_2d_folded); + + RUN_TEST(test_subtract_interior_mean); + RUN_TEST(test_is_pow2_plus1); + + return UNITY_END(); +} From e85a54b0e4929a3e01aa538eb402aeabffdfaac7 Mon Sep 17 00:00:00 2001 From: shaia Date: Sun, 12 Jul 2026 09:27:13 +0300 Subject: [PATCH 2/4] Compile multigrid transfer TU into test only for shared builds The mg_* operators are already part of CFD::Library in static builds (the default here), so compiling multigrid_transfer.c into the test binary risks duplicate definitions; only shared builds, where the symbols are hidden, need the direct compilation. Addresses PR #200 review feedback. --- CMakeLists.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ac04e8c7..b9c9a6b6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -351,9 +351,13 @@ 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 - lib/src/solvers/linear/cpu/multigrid_transfer.c) + $<$:${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 From 47b1a2ac7ef432809bcbd82d7fbb8a540f83f3d0 Mon Sep 17 00:00:00 2001 From: shaia Date: Sun, 12 Jul 2026 09:30:17 +0300 Subject: [PATCH 3/4] Set last-status when poisson_solver_create returns NULL Callers that only see NULL (e.g. the error-handling demo reading cfd_get_last_status) previously got a stale status. Report CFD_ERROR_UNSUPPORTED for unavailable backends and CFD_ERROR_INVALID for unknown methods, matching the factory convention in cfd_solver_create. The non-multigrid branches had the same pre-existing gap and are fixed uniformly. Addresses PR #200 review feedback. --- lib/src/solvers/linear/linear_solver.c | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/lib/src/solvers/linear/linear_solver.c b/lib/src/solvers/linear/linear_solver.c index 8f8d1776..5ba4d0dd 100644 --- a/lib/src/solvers/linear/linear_solver.c +++ b/lib/src/solvers/linear/linear_solver.c @@ -155,6 +155,15 @@ static poisson_solver_backend_t select_best_backend(void) { return POISSON_BACKEND_SCALAR; } +/* Factories must set a specific last-status on NULL (see cfd_solver_create), + * so callers can distinguish an unavailable backend from allocation failure. */ +static poisson_solver_t* backend_unavailable(const char* method_name) { + char msg[128]; + snprintf(msg, sizeof(msg), "Requested backend not available for %s", method_name); + cfd_set_error(CFD_ERROR_UNSUPPORTED, msg); + return NULL; +} + poisson_solver_t* poisson_solver_create( poisson_solver_method_t method, poisson_solver_backend_t backend) @@ -179,7 +188,7 @@ poisson_solver_t* poisson_solver_create( case POISSON_BACKEND_SCALAR: return create_jacobi_scalar_solver(); default: - return NULL; /* Requested backend not available for Jacobi */ + return backend_unavailable("Jacobi"); } case POISSON_METHOD_SOR: @@ -194,7 +203,7 @@ poisson_solver_t* poisson_solver_create( case POISSON_BACKEND_SCALAR: return create_sor_scalar_solver(); default: - return NULL; /* SOR not available for OMP */ + return backend_unavailable("SOR"); } case POISSON_METHOD_REDBLACK_SOR: @@ -212,7 +221,7 @@ poisson_solver_t* poisson_solver_create( case POISSON_BACKEND_SCALAR: return create_redblack_scalar_solver(); default: - return NULL; /* Requested backend not available for Red-Black */ + return backend_unavailable("Red-Black SOR"); } case POISSON_METHOD_CG: @@ -230,7 +239,7 @@ poisson_solver_t* poisson_solver_create( case POISSON_BACKEND_SCALAR: return create_cg_scalar_solver(); default: - return NULL; /* Requested backend not available for CG */ + return backend_unavailable("CG"); } case POISSON_METHOD_BICGSTAB: @@ -244,7 +253,7 @@ poisson_solver_t* poisson_solver_create( case POISSON_BACKEND_SCALAR: return create_bicgstab_scalar_solver(); default: - return NULL; /* Requested backend not available for BiCGSTAB */ + return backend_unavailable("BiCGSTAB"); } case POISSON_METHOD_GMRES: @@ -258,7 +267,7 @@ poisson_solver_t* poisson_solver_create( case POISSON_BACKEND_SCALAR: return create_gmres_scalar_solver(); default: - return NULL; /* Requested backend not available for GMRES */ + return backend_unavailable("GMRES"); } case POISSON_METHOD_MULTIGRID: @@ -269,9 +278,10 @@ poisson_solver_t* poisson_solver_create( backend == POISSON_BACKEND_SCALAR) { return create_multigrid_scalar_solver(); } - return NULL; + return backend_unavailable("multigrid"); default: + cfd_set_error(CFD_ERROR_INVALID, "Unknown Poisson solver method"); return NULL; } } From b9c9c3ab8c4b4f2efa2316cfdd9eaef3f28290ca Mon Sep 17 00:00:00 2001 From: shaia Date: Sun, 12 Jul 2026 09:41:45 +0300 Subject: [PATCH 4/4] Clarify linkage wording for multigrid transfer operators MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The operators are not declared static — they have external linkage and are merely not exported from the shared library. Saying "internal linkage" misstated how the unit test is able to build against them. Addresses PR #200 review feedback. --- lib/src/solvers/linear/multigrid_internal.h | 9 ++++++--- tests/math/test_multigrid_operators.c | 6 ++++-- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/lib/src/solvers/linear/multigrid_internal.h b/lib/src/solvers/linear/multigrid_internal.h index ec769a38..3c841256 100644 --- a/lib/src/solvers/linear/multigrid_internal.h +++ b/lib/src/solvers/linear/multigrid_internal.h @@ -3,8 +3,10 @@ * @brief Internal data structures and grid-transfer operators for geometric multigrid * * Not part of the public API. The transfer operators live in - * cpu/multigrid_transfer.c with internal (non-exported) linkage so the - * operator unit test can compile that translation unit directly. + * cpu/multigrid_transfer.c with external linkage but are not exported from + * the shared library (no CFD_LIBRARY_EXPORT, hidden visibility): the operator + * unit test links them from the static library, or compiles the translation + * unit directly in shared builds. */ #ifndef CFD_MULTIGRID_INTERNAL_H @@ -77,7 +79,8 @@ static inline bool mg_is_pow2_plus1(size_t n) { /* ============================================================================ * GRID-TRANSFER OPERATORS (cpu/multigrid_transfer.c) * - * Internal linkage across the library; intentionally NOT CFD_LIBRARY_EXPORT. + * External linkage, intentionally NOT CFD_LIBRARY_EXPORT (hidden outside a + * shared library build). * All write coarse/fine INTERIOR points only. Restriction never reads fine * boundary values (coarse-interior stencils stay within the fine interior); * prolongation reads coarse boundary values, which the caller must have set diff --git a/tests/math/test_multigrid_operators.c b/tests/math/test_multigrid_operators.c index 5114d2d3..0a476bb9 100644 --- a/tests/math/test_multigrid_operators.c +++ b/tests/math/test_multigrid_operators.c @@ -2,8 +2,10 @@ * @file test_multigrid_operators.c * @brief Unit tests for multigrid grid-transfer operators * - * Compiles lib/src/solvers/linear/cpu/multigrid_transfer.c directly (the - * operators have internal linkage and are not exported from the library). + * Exercises lib/src/solvers/linear/cpu/multigrid_transfer.c. The operators + * have external linkage but are not exported from the shared library, so + * static builds link them from CFD::Library while shared builds compile the + * translation unit into this test (see CMakeLists.txt). * Deliberately uses no poisson_solver_* API and no cfd runtime — pure math * verification of the operators before they are used in the solver. *