From ad2a8a39574848d00a51dde7f0c232139a28882a Mon Sep 17 00:00:00 2001 From: shaia Date: Fri, 17 Jul 2026 17:48:22 +0300 Subject: [PATCH 1/9] Add multigrid-preconditioned CG to the scalar Poisson solver One symmetric V(2,2) weighted-Jacobi multigrid cycle in Dirichlet mode per preconditioner apply: Dirichlet matches the interior-only operator CG applies to its zero-halo Krylov vectors (keeping M nonsingular), and equal pre/post Jacobi sweeps keep M symmetric as CG requires. Gives grid-size-independent convergence (5 iterations at 33..129 vs 50-170 for plain CG at tol 1e-8) on uniform grids where Jacobi preconditioning gives no benefit. The convenience API gains a POISSON_SOLVER_PCG_MG_SCALAR preset because poisson_solve_3d has no params argument: the preset carries the preconditioner into its per-preset cached instance. Backends without the preconditioner (OMP/SIMD CG, GMRES, GPU CG) reject POISSON_PRECOND_MULTIGRID explicitly instead of silently ignoring it. --- CMakeLists.txt | 3 + docs/guides/examples.md | 13 +- examples/poisson_solver_tuning.c | 9 +- lib/include/cfd/solvers/poisson_solver.h | 10 +- .../linear/avx2/linear_solver_cg_avx2.c | 5 + lib/src/solvers/linear/cpu/linear_solver_cg.c | 110 +++- .../solvers/linear/cpu/linear_solver_gmres.c | 5 + .../linear/gpu/poisson_solver_cg_gpu.cu | 6 +- lib/src/solvers/linear/linear_solver.c | 23 +- .../solvers/linear/linear_solver_internal.h | 15 + .../linear/neon/linear_solver_cg_neon.c | 5 + .../solvers/linear/omp/linear_solver_cg_omp.c | 5 + .../linear/omp/linear_solver_gmres_omp.c | 5 + .../linear_solver_gmres_simd_template.h | 5 + tests/math/test_mg_pcg_convergence.c | 512 ++++++++++++++++++ 15 files changed, 714 insertions(+), 17 deletions(-) create mode 100644 tests/math/test_mg_pcg_convergence.c diff --git a/CMakeLists.txt b/CMakeLists.txt index b9c9a6b6..6d161b7f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -348,6 +348,7 @@ if(BUILD_TESTS) add_executable(test_divergence_free tests/math/test_divergence_free.c) add_executable(test_mms tests/math/test_mms.c) add_executable(test_pcg_convergence tests/math/test_pcg_convergence.c) + add_executable(test_mg_pcg_convergence tests/math/test_mg_pcg_convergence.c) add_executable(test_residual_computation tests/math/test_residual_computation.c) add_executable(test_cg_scaling tests/math/test_cg_scaling.c) add_executable(test_nonuniform_grid tests/math/test_nonuniform_grid.c) @@ -478,6 +479,7 @@ if(BUILD_TESTS) target_link_libraries(test_divergence_free PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_mms PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_pcg_convergence PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_mg_pcg_convergence PRIVATE CFD::Library unity $<$>:m>) 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>) @@ -642,6 +644,7 @@ if(BUILD_TESTS) add_test(NAME DivergenceFreeTest COMMAND test_divergence_free) add_test(NAME MmsTest COMMAND test_mms) add_test(NAME PCGConvergenceTest COMMAND test_pcg_convergence) + add_test(NAME MgPcgConvergenceTest COMMAND test_mg_pcg_convergence) add_test(NAME ResidualComputationTest COMMAND test_residual_computation) add_test(NAME CgScalingTest COMMAND test_cg_scaling) add_test(NAME NonuniformGridTest COMMAND test_nonuniform_grid) diff --git a/docs/guides/examples.md b/docs/guides/examples.md index 413b6491..2fa2584b 100644 --- a/docs/guides/examples.md +++ b/docs/guides/examples.md @@ -681,12 +681,12 @@ CFD Platform Diagnostics - Error handling for unavailable solvers **Sections:** -1. Method comparison (Jacobi, SOR, Red-Black SOR, CG, CG+Jacobi PC, BiCGSTAB) on scalar backend +1. Method comparison (Jacobi, SOR, Red-Black SOR, CG, CG+Jacobi PC, CG+Multigrid PC, BiCGSTAB) on scalar backend 2. Backend comparison (CG on Scalar, SIMD, OMP) 3. Convenience API demo (`poisson_solve()`) 4. Error handling (requesting multigrid on an unavailable backend — GPU) -**Problem:** Solves ∇²p = -2π²sin(πx)sin(πy) on a 64×64 grid using the library's default homogeneous Neumann boundary conditions. The reported L2 error compares methods/backends against a common reference field — not against the Dirichlet analytical solution sin(πx)sin(πy), since BCs differ. +**Problem:** Solves ∇²p = -2π²sin(πx)sin(πy) on a 65×65 grid (2^k+1, so the multigrid preconditioner can build its hierarchy) using the library's default homogeneous Neumann boundary conditions. The reported L2 error compares methods/backends against a common reference field — not against the Dirichlet analytical solution sin(πx)sin(πy), since BCs differ. Standalone multigrid is omitted: this RHS has a nonzero interior mean, which the true Neumann system cannot converge on (the same reason the stationary methods report max_iter). **Run:** ```bash @@ -697,10 +697,11 @@ CFD Platform Diagnostics ``` --- Method Comparison (Scalar Backend) --- Method Iters Residual L2 Error Time Status - Jacobi 10001 res=8.3e+00 L2=4.8e+00 5548 ms max_iter - CG 1 res=5.1e-11 L2=1.1e-04 2 ms converged - CG + Jacobi PC 1 res=5.9e-11 L2=1.1e-04 2 ms converged - BiCGSTAB 1 res=5.1e-11 L2=1.1e-04 2 ms converged + Jacobi 10000 res=8.3e+00 L2=-1.0e+00 394 ms max_iter + CG 1 res=6.3e-11 L2=1.0e-04 0 ms converged + CG + Jacobi PC 1 res=6.3e-11 L2=1.0e-04 0 ms converged + CG + Multigrid PC 7 res=1.4e-06 L2=1.0e-04 2 ms converged + BiCGSTAB 1 res=6.3e-11 L2=1.0e-04 0 ms converged ``` --- diff --git a/examples/poisson_solver_tuning.c b/examples/poisson_solver_tuning.c index 5d507871..95879720 100644 --- a/examples/poisson_solver_tuning.c +++ b/examples/poisson_solver_tuning.c @@ -124,7 +124,8 @@ int main(void) { printf("Poisson Solver Tuning Example\n"); printf("=============================\n"); - size_t nx = 64, ny = 64; + /* 65 = 2^6+1 so the multigrid rows below can build their hierarchy */ + size_t nx = 65, ny = 65; double dx = 1.0 / (double)(nx - 1); double dy = 1.0 / (double)(ny - 1); @@ -163,8 +164,14 @@ int main(void) { POISSON_PRECOND_NONE, nx, ny, dx, dy, rhs, p, p_temp, p_exact); benchmark_method("CG + Jacobi PC", POISSON_METHOD_CG, POISSON_BACKEND_SCALAR, POISSON_PRECOND_JACOBI, nx, ny, dx, dy, rhs, p, p_temp, p_exact); + benchmark_method("CG + Multigrid PC", POISSON_METHOD_CG, POISSON_BACKEND_SCALAR, + POISSON_PRECOND_MULTIGRID, nx, ny, dx, dy, rhs, p, p_temp, p_exact); benchmark_method("BiCGSTAB", POISSON_METHOD_BICGSTAB, POISSON_BACKEND_SCALAR, POISSON_PRECOND_NONE, nx, ny, dx, dy, rhs, p, p_temp, p_exact); + /* Standalone POISSON_METHOD_MULTIGRID is omitted here: this example's RHS + * has a nonzero interior mean, which the true Neumann system cannot + * converge on (same reason the stationary methods above report max_iter). + * See tests/math/test_multigrid_convergence.c for a compatible setup. */ /* Section 2: Backend comparison (CG method) */ printf("\n--- Backend Comparison (CG Method) ---\n"); diff --git a/lib/include/cfd/solvers/poisson_solver.h b/lib/include/cfd/solvers/poisson_solver.h index 40d892e2..1887ab62 100644 --- a/lib/include/cfd/solvers/poisson_solver.h +++ b/lib/include/cfd/solvers/poisson_solver.h @@ -87,8 +87,10 @@ typedef enum { * Preconditioner types for iterative solvers */ typedef enum { - POISSON_PRECOND_NONE = 0, /**< No preconditioning (default) */ - POISSON_PRECOND_JACOBI = 1 /**< Diagonal (Jacobi) preconditioning */ + POISSON_PRECOND_NONE = 0, /**< No preconditioning (default) */ + POISSON_PRECOND_JACOBI = 1, /**< Diagonal (Jacobi) preconditioning */ + POISSON_PRECOND_MULTIGRID = 2 /**< One geometric-multigrid V-cycle per apply + (scalar CG only; grid dims must be 2^k+1) */ } poisson_precond_type_t; /** @@ -461,7 +463,9 @@ typedef enum { 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_MG_SCALAR = 9 /**< Geometric multigrid with scalar backend (grid dims must be 2^k+1) */ + POISSON_SOLVER_MG_SCALAR = 9, /**< Geometric multigrid with scalar backend (grid dims must be 2^k+1) */ + POISSON_SOLVER_PCG_MG_SCALAR = 10 /**< CG with multigrid V-cycle preconditioner, 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/avx2/linear_solver_cg_avx2.c b/lib/src/solvers/linear/avx2/linear_solver_cg_avx2.c index a3a141df..8664b014 100644 --- a/lib/src/solvers/linear/avx2/linear_solver_cg_avx2.c +++ b/lib/src/solvers/linear/avx2/linear_solver_cg_avx2.c @@ -393,6 +393,11 @@ static cfd_status_t cg_avx2_init( double dx, double dy, double dz, const poisson_solver_params_t* params) { + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + /* Use aligned allocation for SIMD context */ cg_avx2_context_t* ctx = (cg_avx2_context_t*)cfd_aligned_calloc( 1, sizeof(cg_avx2_context_t)); diff --git a/lib/src/solvers/linear/cpu/linear_solver_cg.c b/lib/src/solvers/linear/cpu/linear_solver_cg.c index 0afc0033..08bd7854 100644 --- a/lib/src/solvers/linear/cpu/linear_solver_cg.c +++ b/lib/src/solvers/linear/cpu/linear_solver_cg.c @@ -53,6 +53,9 @@ typedef struct { double* p; /* Search direction */ double* Ap; /* A * p (Laplacian applied to p) */ + poisson_solver_t* mg_precond; /* Inner MG V-cycle solver (NULL unless + POISSON_PRECOND_MULTIGRID) */ + poisson_precond_type_t precond_type; /* Which preconditioner use_precond refers to */ int use_precond; /* Flag: is preconditioner enabled? */ int initialized; } cg_context_t; @@ -192,6 +195,38 @@ static void apply_jacobi_precond(const double* r, double* z, } } +/** + * Apply multigrid preconditioner: z = M^{-1} * r for A = -nabla^2. + * + * The inner solver runs one V-cycle on nabla^2 z' = r from a zero guess; + * by linearity of the fixed cycle, z = -z' approximates (-nabla^2)^{-1} r. + * r can be passed directly as the MG rhs: its boundary stays zero and the + * MG smoothers/residual never read rhs boundary values. The cycle uses + * Dirichlet mode with symmetric Jacobi smoothing (set at init), matching + * the zero halo that CG's interior-only updates impose on Krylov vectors + * and keeping M symmetric positive definite. + */ +static cfd_status_t apply_mg_precond(poisson_solver_t* mg, const double* r, double* z, + size_t n_total, size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + memset(z, 0, n_total * sizeof(double)); + + cfd_status_t status = poisson_solver_iterate(mg, z, NULL, r, NULL); + if (status != CFD_SUCCESS) { + return status; + } + + 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++) { + size_t idx = (k * stride_z) + (IDX_2D(i, j, nx)); + z[idx] = -z[idx]; + } + } + } + return CFD_SUCCESS; +} + /* ============================================================================ * CG SCALAR IMPLEMENTATION * ============================================================================ */ @@ -216,7 +251,9 @@ static cfd_status_t cg_scalar_init( ctx->diag_inv = 1.0 / (2.0 / ctx->dx2 + 2.0 / ctx->dy2 + 2.0 * ctx->inv_dz2); /* Check if preconditioner is enabled */ - ctx->use_precond = (params && params->preconditioner == POISSON_PRECOND_JACOBI); + ctx->precond_type = params ? params->preconditioner : POISSON_PRECOND_NONE; + ctx->use_precond = (ctx->precond_type == POISSON_PRECOND_JACOBI || + ctx->precond_type == POISSON_PRECOND_MULTIGRID); /* Allocate working vectors */ size_t n = nx * ny * nz; @@ -245,6 +282,39 @@ static cfd_status_t cg_scalar_init( } } + if (ctx->precond_type == POISSON_PRECOND_MULTIGRID) { + ctx->mg_precond = create_multigrid_scalar_solver(); + if (!ctx->mg_precond) { + cfd_free(ctx->r); + cfd_free(ctx->z); + cfd_free(ctx->p); + cfd_free(ctx->Ap); + cfd_free(ctx); + return CFD_ERROR_NOMEM; + } + + /* One V-cycle per apply. Dirichlet mode matches CG's interior + * operator (Krylov vectors carry a permanent zero halo) and keeps M + * nonsingular; Jacobi smoothing with equal pre/post sweeps keeps M + * symmetric, as CG requires. */ + poisson_solver_params_t mg_params = poisson_solver_params_default(); + mg_params.mg_cycle = MG_CYCLE_V; + mg_params.mg_smoother = MG_SMOOTHER_JACOBI; + mg_params.mg_bc = MG_BC_DIRICHLET; + + cfd_status_t mg_status = poisson_solver_init( + ctx->mg_precond, nx, ny, nz, dx, dy, dz, &mg_params); + if (mg_status != CFD_SUCCESS) { + poisson_solver_destroy(ctx->mg_precond); + cfd_free(ctx->r); + cfd_free(ctx->z); + cfd_free(ctx->p); + cfd_free(ctx->Ap); + cfd_free(ctx); + return mg_status; /* CFD_ERROR_INVALID for non-2^k+1 dims */ + } + } + ctx->initialized = 1; solver->context = ctx; return CFD_SUCCESS; @@ -253,6 +323,7 @@ static cfd_status_t cg_scalar_init( static void cg_scalar_destroy(poisson_solver_t* solver) { if (solver && solver->context) { cg_context_t* ctx = (cg_context_t*)solver->context; + poisson_solver_destroy(ctx->mg_precond); cfd_free(ctx->r); cfd_free(ctx->z); cfd_free(ctx->p); @@ -322,11 +393,28 @@ static cfd_status_t cg_scalar_solve( /* Compute initial residual: r_0 = b - A*x_0 */ compute_residual(x, rhs, r, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + size_t n_total = nx * ny * solver->nz; + double initial_res = sqrt(dot_product(r, r, nx, ny, k_start, k_end, stride_z)); + /* Initialize search direction and rho */ double rho; if (use_precond) { /* z_0 = M^{-1} r_0 */ - apply_jacobi_precond(r, z, nx, ny, diag_inv, k_start, k_end, stride_z); + if (ctx->precond_type == POISSON_PRECOND_MULTIGRID) { + cfd_status_t precond_status = apply_mg_precond( + ctx->mg_precond, r, z, n_total, nx, ny, k_start, k_end, stride_z); + if (precond_status != CFD_SUCCESS) { + if (stats) { + stats->status = POISSON_ERROR; + stats->iterations = 0; + stats->final_residual = initial_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return precond_status; + } + } else { + apply_jacobi_precond(r, z, nx, ny, diag_inv, k_start, k_end, stride_z); + } /* p_0 = z_0 */ copy_vector(z, p, nx, ny, k_start, k_end, stride_z); /* rho_0 = (r_0, z_0) */ @@ -338,8 +426,6 @@ static cfd_status_t cg_scalar_solve( rho = dot_product(r, r, nx, ny, k_start, k_end, stride_z); } - double initial_res = sqrt(dot_product(r, r, nx, ny, k_start, k_end, stride_z)); - if (stats) { stats->initial_residual = initial_res; } @@ -386,7 +472,21 @@ static cfd_status_t cg_scalar_solve( double rho_new; if (use_precond) { /* z_{k+1} = M^{-1} r_{k+1} */ - apply_jacobi_precond(r, z, nx, ny, diag_inv, k_start, k_end, stride_z); + if (ctx->precond_type == POISSON_PRECOND_MULTIGRID) { + cfd_status_t precond_status = apply_mg_precond( + ctx->mg_precond, r, z, n_total, nx, ny, k_start, k_end, stride_z); + if (precond_status != CFD_SUCCESS) { + if (stats) { + stats->status = POISSON_ERROR; + stats->iterations = iter + 1; + stats->final_residual = res_norm; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return precond_status; + } + } else { + apply_jacobi_precond(r, z, nx, ny, diag_inv, k_start, k_end, stride_z); + } /* rho_new = (r_{k+1}, z_{k+1}) */ rho_new = dot_product(r, z, nx, ny, k_start, k_end, stride_z); } else { diff --git a/lib/src/solvers/linear/cpu/linear_solver_gmres.c b/lib/src/solvers/linear/cpu/linear_solver_gmres.c index 650cafe2..6baeb21f 100644 --- a/lib/src/solvers/linear/cpu/linear_solver_gmres.c +++ b/lib/src/solvers/linear/cpu/linear_solver_gmres.c @@ -314,6 +314,11 @@ static cfd_status_t gmres_scalar_init( double dx, double dy, double dz, const poisson_solver_params_t* params) { + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + gmres_context_t* ctx = (gmres_context_t*)cfd_calloc(1, sizeof(gmres_context_t)); if (!ctx) { return CFD_ERROR_NOMEM; diff --git a/lib/src/solvers/linear/gpu/poisson_solver_cg_gpu.cu b/lib/src/solvers/linear/gpu/poisson_solver_cg_gpu.cu index b07ddbbd..6781ecf8 100644 --- a/lib/src/solvers/linear/gpu/poisson_solver_cg_gpu.cu +++ b/lib/src/solvers/linear/gpu/poisson_solver_cg_gpu.cu @@ -61,7 +61,11 @@ static cfd_status_t cg_gpu_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) { - (void)params; + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + if (!gpu_is_available()) { cfd_set_error(CFD_ERROR_UNSUPPORTED, "CUDA GPU not available at runtime"); return CFD_ERROR_UNSUPPORTED; diff --git a/lib/src/solvers/linear/linear_solver.c b/lib/src/solvers/linear/linear_solver.c index 5ba4d0dd..1c845743 100644 --- a/lib/src/solvers/linear/linear_solver.c +++ b/lib/src/solvers/linear/linear_solver.c @@ -588,6 +588,7 @@ 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; +static poisson_solver_t* g_cached_pcg_mg_scalar = NULL; /** * Cleanup cached solvers (called at program exit) @@ -633,6 +634,10 @@ static void cleanup_cached_solvers(void) { poisson_solver_destroy(g_cached_mg_scalar); g_cached_mg_scalar = NULL; } + if (g_cached_pcg_mg_scalar) { + poisson_solver_destroy(g_cached_pcg_mg_scalar); + g_cached_pcg_mg_scalar = NULL; + } } int poisson_solve_3d( @@ -706,6 +711,12 @@ int poisson_solve_3d( backend = POISSON_BACKEND_SCALAR; break; + case POISSON_SOLVER_PCG_MG_SCALAR: + solver_ptr = &g_cached_pcg_mg_scalar; + method = POISSON_METHOD_CG; + backend = POISSON_BACKEND_SCALAR; + break; + default: CFD_LOG_ERROR("poisson", "poisson_solve_3d: Unknown solver type %d", solver_type); return -1; @@ -737,8 +748,18 @@ int poisson_solve_3d( *solver_ptr = poisson_solver_create(method, backend); if (*solver_ptr) { + /* The convenience API has no params argument, so the PCG_MG + * preset carries its preconditioner into the cached instance. */ + poisson_solver_params_t pcg_mg_params; + const poisson_solver_params_t* init_params = NULL; + if (solver_type == POISSON_SOLVER_PCG_MG_SCALAR) { + pcg_mg_params = poisson_solver_params_default(); + pcg_mg_params.preconditioner = POISSON_PRECOND_MULTIGRID; + init_params = &pcg_mg_params; + } + cfd_status_t init_status = - poisson_solver_init(*solver_ptr, nx, ny, nz, dx, dy, dz, NULL); + poisson_solver_init(*solver_ptr, nx, ny, nz, dx, dy, dz, init_params); 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 fa4b3f9c..9d003f7d 100644 --- a/lib/src/solvers/linear/linear_solver_internal.h +++ b/lib/src/solvers/linear/linear_solver_internal.h @@ -73,6 +73,21 @@ poisson_solver_t* create_gmres_omp_solver(void); /* Geometric multigrid solver (scalar backend only; V/W/F cycles) */ poisson_solver_t* create_multigrid_scalar_solver(void); +/** + * Reject POISSON_PRECOND_MULTIGRID on backends that don't implement it. + * Only the scalar CG solver supports the MG preconditioner; silently + * ignoring it would be a forbidden silent fallback. + */ +static inline cfd_status_t poisson_solver_reject_mg_precond( + const poisson_solver_params_t* params) { + if (params && params->preconditioner == POISSON_PRECOND_MULTIGRID) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, + "POISSON_PRECOND_MULTIGRID is only supported by the scalar CG solver"); + return CFD_ERROR_UNSUPPORTED; + } + return CFD_SUCCESS; +} + /* ============================================================================ * CG ALGORITHM CONSTANTS * ============================================================================ */ diff --git a/lib/src/solvers/linear/neon/linear_solver_cg_neon.c b/lib/src/solvers/linear/neon/linear_solver_cg_neon.c index b04543cd..ca90f7f3 100644 --- a/lib/src/solvers/linear/neon/linear_solver_cg_neon.c +++ b/lib/src/solvers/linear/neon/linear_solver_cg_neon.c @@ -389,6 +389,11 @@ static cfd_status_t cg_neon_init( double dx, double dy, double dz, const poisson_solver_params_t* params) { + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + cg_neon_context_t* ctx = (cg_neon_context_t*)cfd_aligned_calloc( 1, sizeof(cg_neon_context_t)); if (!ctx) { diff --git a/lib/src/solvers/linear/omp/linear_solver_cg_omp.c b/lib/src/solvers/linear/omp/linear_solver_cg_omp.c index 3b340afd..fbd24c48 100644 --- a/lib/src/solvers/linear/omp/linear_solver_cg_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_cg_omp.c @@ -203,6 +203,11 @@ static cfd_status_t cg_omp_init( double dx, double dy, double dz, const poisson_solver_params_t* params) { + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + cg_omp_context_t* ctx = (cg_omp_context_t*)cfd_calloc(1, sizeof(cg_omp_context_t)); if (!ctx) { return CFD_ERROR_NOMEM; diff --git a/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c b/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c index 3937c7e0..26fcecc9 100644 --- a/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_gmres_omp.c @@ -282,6 +282,11 @@ static cfd_status_t gmres_omp_init( double dx, double dy, double dz, const poisson_solver_params_t* params) { + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + gmres_omp_context_t* ctx = (gmres_omp_context_t*)cfd_calloc(1, sizeof(gmres_omp_context_t)); if (!ctx) { return CFD_ERROR_NOMEM; diff --git a/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h b/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h index 4266119d..e4f1ce2b 100644 --- a/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h +++ b/lib/src/solvers/linear/simd_template/linear_solver_gmres_simd_template.h @@ -411,6 +411,11 @@ static cfd_status_t SIMD_FUNC(gmres_init)( double dx, double dy, double dz, const poisson_solver_params_t* params) { + cfd_status_t precond_status = poisson_solver_reject_mg_precond(params); + if (precond_status != CFD_SUCCESS) { + return precond_status; + } + gmres_simd_context_t* ctx = (gmres_simd_context_t*)cfd_aligned_calloc( 1, sizeof(gmres_simd_context_t)); if (!ctx) { diff --git a/tests/math/test_mg_pcg_convergence.c b/tests/math/test_mg_pcg_convergence.c new file mode 100644 index 00000000..bdad60f2 --- /dev/null +++ b/tests/math/test_mg_pcg_convergence.c @@ -0,0 +1,512 @@ +/** + * @file test_mg_pcg_convergence.c + * @brief Multigrid-preconditioned Conjugate Gradient convergence tests + * + * These tests verify that CG with the multigrid V-cycle preconditioner + * (POISSON_PRECOND_MULTIGRID): + * - Converges to the same solution as standard CG + * - Needs a grid-size-independent iteration count (the MG payoff), + * far fewer iterations than plain CG on larger grids + * - Works in 3D + * - Rejects non-2^k+1 grid dimensions at init + * - Is explicitly rejected (CFD_ERROR_UNSUPPORTED) by non-scalar CG and + * by GMRES, instead of being silently ignored + * - Is reachable through the poisson_solve_3d() convenience preset + * POISSON_SOLVER_PCG_MG_SCALAR with the legacy return contract + * + * The preconditioner runs one symmetric V(2,2) weighted-Jacobi cycle in + * Dirichlet mode per apply (see linear_solver_cg.c), so M is SPD and CG + * theory applies. + */ + +#include "unity.h" +#include "cfd/solvers/poisson_solver.h" +#include "cfd/core/memory.h" +#include "cfd/core/indexing.h" +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +/* ============================================================================ + * TEST PARAMETERS + * ============================================================================ */ + +#define DOMAIN_MIN 0.0 +#define DOMAIN_MAX 1.0 + +#define TOLERANCE 1e-8 +#define MAX_ITERATIONS 5000 + +/* Grid-independence bound for MG-PCG: a working MG preconditioner keeps the + * outer CG count roughly constant across grid sizes. */ +#define MG_PCG_MAX_ITERS 25 + +/* On larger grids plain CG needs O(n) iterations while MG-PCG stays flat; + * require at least a 2x reduction from 65x65 up. */ +#define MG_PCG_MIN_SPEEDUP 2.0 + +/* ============================================================================ + * HELPER FUNCTIONS + * ============================================================================ */ + +static double* create_field_3d(size_t nx, size_t ny, size_t nz) { + return (double*)cfd_calloc(nx * ny * nz, sizeof(double)); +} + +static void init_nontrivial_guess(double* p, size_t nx, size_t ny) { + for (size_t j = 0; j < ny; j++) { + for (size_t i = 0; i < nx; i++) { + p[IDX_2D(i, j, nx)] = ((i + j) % 2 == 0) ? 1.0 : -1.0; + } + } +} + +/** + * Initialize sinusoidal RHS compatible with Neumann BCs. + * f(x,y) = cos(2πx)cos(2πy) with discrete interior mean subtracted. + */ +static void init_sinusoidal_rhs(double* rhs, size_t nx, size_t ny, + double dx, double dy) { + for (size_t j = 0; j < ny; j++) { + double y = DOMAIN_MIN + j * dy; + for (size_t i = 0; i < nx; i++) { + double x = DOMAIN_MIN + i * dx; + rhs[IDX_2D(i, j, nx)] = cos(2.0 * M_PI * x) * cos(2.0 * M_PI * y); + } + } + + double interior_sum = 0.0; + size_t interior_count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + interior_sum += rhs[IDX_2D(i, j, nx)]; + interior_count++; + } + } + + if (interior_count > 0) { + double interior_mean = interior_sum / (double)interior_count; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + rhs[IDX_2D(i, j, nx)] -= interior_mean; + } + } + } + + for (size_t i = 0; i < nx; i++) { + rhs[i] = 0.0; + rhs[(ny - 1) * nx + i] = 0.0; + } + for (size_t j = 1; j < ny - 1; j++) { + rhs[j * nx] = 0.0; + rhs[j * nx + (nx - 1)] = 0.0; + } +} + +/** + * 3D analog: f = cos(2πx)cos(2πy)cos(2πz), interior mean subtracted, + * boundary RHS zeroed. + */ +static void init_sinusoidal_rhs_3d(double* rhs, size_t nx, size_t ny, size_t nz, + double dx, double dy, double dz) { + size_t plane = nx * ny; + for (size_t k = 0; k < nz; k++) { + double z = DOMAIN_MIN + k * dz; + for (size_t j = 0; j < ny; j++) { + double y = DOMAIN_MIN + j * dy; + for (size_t i = 0; i < nx; i++) { + double x = DOMAIN_MIN + i * dx; + rhs[k * plane + IDX_2D(i, j, nx)] = + cos(2.0 * M_PI * x) * cos(2.0 * M_PI * y) * cos(2.0 * M_PI * z); + } + } + } + + double interior_sum = 0.0; + size_t interior_count = 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++) { + interior_sum += rhs[k * plane + IDX_2D(i, j, nx)]; + interior_count++; + } + } + } + double interior_mean = (interior_count > 0) + ? interior_sum / (double)interior_count : 0.0; + + for (size_t k = 0; k < nz; k++) { + for (size_t j = 0; j < ny; j++) { + for (size_t i = 0; i < nx; i++) { + size_t idx = k * plane + IDX_2D(i, j, nx); + int interior = (k > 0 && k < nz - 1 && j > 0 && j < ny - 1 && + i > 0 && i < nx - 1); + rhs[idx] = interior ? rhs[idx] - interior_mean : 0.0; + } + } + } +} + +/** L2 norm of the interior difference of two fields (3D-aware; nz=1 for 2D) */ +static double compute_l2_difference_3d(const double* a, const double* b, + size_t nx, size_t ny, size_t nz) { + size_t plane = nx * ny; + size_t k_start = (nz > 1) ? 1 : 0; + size_t k_end = (nz > 1) ? nz - 1 : 1; + double sum_sq = 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++) { + size_t idx = k * plane + IDX_2D(i, j, nx); + double diff = a[idx] - b[idx]; + sum_sq += diff * diff; + count++; + } + } + } + return sqrt(sum_sq / (double)count); +} + +/** + * Subtract the interior mean (Neumann solutions are defined up to a constant; + * CG and MG-PCG may settle on different constants). + */ +static void subtract_interior_mean_3d(double* f, size_t nx, size_t ny, size_t nz) { + size_t plane = nx * ny; + 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 * plane + IDX_2D(i, j, nx)]; + count++; + } + } + } + double mean = sum / (double)count; + 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 * plane + IDX_2D(i, j, nx)] -= mean; + } + } + } +} + +/** + * Run one CG solve on the given problem with the given preconditioner. + * Returns iteration count; asserts convergence. + */ +static int run_cg_solve(double* p, double* p_temp, const double* rhs, + size_t nx, size_t ny, size_t nz, + double dx, double dy, double dz, + poisson_precond_type_t precond) { + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_CG, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL_MESSAGE(solver, "Could not create CG solver"); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = MAX_ITERATIONS; + params.preconditioner = precond; + + cfd_status_t status = poisson_solver_init(solver, nx, ny, nz, dx, dy, dz, ¶ms); + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, status); + + poisson_solver_stats_t stats = poisson_solver_stats_default(); + status = poisson_solver_solve(solver, p, p_temp, rhs, &stats); + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL_INT(POISSON_CONVERGED, stats.status); + + poisson_solver_destroy(solver); + return stats.iterations; +} + +/* ============================================================================ + * TEST: MG-PCG CONVERGES TO THE SAME SOLUTION AS CG + * ============================================================================ */ + +void test_mg_pcg_converges_correctly(void) { + printf("\n Testing MG-PCG converges to correct solution...\n"); + + size_t n = 33; + double h = (DOMAIN_MAX - DOMAIN_MIN) / (n - 1); + + double* p_cg = create_field_3d(n, n, 1); + double* p_pcg = create_field_3d(n, n, 1); + double* p_temp = create_field_3d(n, n, 1); + double* rhs = create_field_3d(n, n, 1); + TEST_ASSERT_NOT_NULL_MESSAGE(p_cg, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_pcg, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_temp, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(rhs, "Memory allocation failed"); + + init_sinusoidal_rhs(rhs, n, n, h, h); + + init_nontrivial_guess(p_cg, n, n); + int cg_iters = run_cg_solve(p_cg, p_temp, rhs, n, n, 1, h, h, 0.0, + POISSON_PRECOND_NONE); + + init_nontrivial_guess(p_pcg, n, n); + int pcg_iters = run_cg_solve(p_pcg, p_temp, rhs, n, n, 1, h, h, 0.0, + POISSON_PRECOND_MULTIGRID); + + printf(" CG: %d iterations\n", cg_iters); + printf(" MG-PCG: %d iterations\n", pcg_iters); + + /* Neumann solutions are defined up to a constant */ + subtract_interior_mean_3d(p_cg, n, n, 1); + subtract_interior_mean_3d(p_pcg, n, n, 1); + + double l2_diff = compute_l2_difference_3d(p_cg, p_pcg, n, n, 1); + printf(" L2 difference between CG and MG-PCG solutions: %.2e\n", l2_diff); + + TEST_ASSERT_TRUE_MESSAGE(l2_diff < 1e-6, + "MG-PCG and CG should converge to the same solution"); + + cfd_free(p_cg); + cfd_free(p_pcg); + cfd_free(p_temp); + cfd_free(rhs); +} + +/* ============================================================================ + * TEST: MG-PCG ITERATION COUNT IS GRID-SIZE-INDEPENDENT + * ============================================================================ */ + +void test_mg_pcg_iteration_reduction(void) { + printf("\n Testing MG-PCG grid-independent iteration count...\n"); + + size_t sizes[] = {33, 65, 129}; + int num_sizes = sizeof(sizes) / sizeof(sizes[0]); + + for (int s = 0; s < num_sizes; s++) { + size_t n = sizes[s]; + double h = (DOMAIN_MAX - DOMAIN_MIN) / (n - 1); + + double* p = create_field_3d(n, n, 1); + double* p_temp = create_field_3d(n, n, 1); + double* rhs = create_field_3d(n, n, 1); + TEST_ASSERT_NOT_NULL_MESSAGE(p, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_temp, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(rhs, "Memory allocation failed"); + + init_sinusoidal_rhs(rhs, n, n, h, h); + + init_nontrivial_guess(p, n, n); + int cg_iters = run_cg_solve(p, p_temp, rhs, n, n, 1, h, h, 0.0, + POISSON_PRECOND_NONE); + + init_nontrivial_guess(p, n, n); + int pcg_iters = run_cg_solve(p, p_temp, rhs, n, n, 1, h, h, 0.0, + POISSON_PRECOND_MULTIGRID); + + double speedup = (double)cg_iters / (double)pcg_iters; + printf(" %3zux%-3zu: CG=%4d iters, MG-PCG=%3d iters (%.1fx)\n", + n, n, cg_iters, pcg_iters, speedup); + + TEST_ASSERT_TRUE_MESSAGE(pcg_iters <= MG_PCG_MAX_ITERS, + "MG-PCG iteration count should be grid-size-independent"); + + if (n >= 65) { + TEST_ASSERT_TRUE_MESSAGE(speedup >= MG_PCG_MIN_SPEEDUP, + "MG-PCG should need far fewer iterations than CG on larger grids"); + } + + cfd_free(p); + cfd_free(p_temp); + cfd_free(rhs); + } +} + +/* ============================================================================ + * TEST: MG-PCG IN 3D + * ============================================================================ */ + +void test_mg_pcg_3d(void) { + printf("\n Testing MG-PCG on a 3D grid...\n"); + + size_t n = 17; + double h = (DOMAIN_MAX - DOMAIN_MIN) / (n - 1); + + double* p_cg = create_field_3d(n, n, n); + double* p_pcg = create_field_3d(n, n, n); + double* p_temp = create_field_3d(n, n, n); + double* rhs = create_field_3d(n, n, n); + TEST_ASSERT_NOT_NULL_MESSAGE(p_cg, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_pcg, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_temp, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(rhs, "Memory allocation failed"); + + init_sinusoidal_rhs_3d(rhs, n, n, n, h, h, h); + + int cg_iters = run_cg_solve(p_cg, p_temp, rhs, n, n, n, h, h, h, + POISSON_PRECOND_NONE); + int pcg_iters = run_cg_solve(p_pcg, p_temp, rhs, n, n, n, h, h, h, + POISSON_PRECOND_MULTIGRID); + + printf(" CG: %d iterations\n", cg_iters); + printf(" MG-PCG: %d iterations\n", pcg_iters); + + subtract_interior_mean_3d(p_cg, n, n, n); + subtract_interior_mean_3d(p_pcg, n, n, n); + + double l2_diff = compute_l2_difference_3d(p_cg, p_pcg, n, n, n); + printf(" L2 difference: %.2e\n", l2_diff); + + TEST_ASSERT_TRUE_MESSAGE(l2_diff < 1e-6, + "3D MG-PCG and CG should converge to the same solution"); + + cfd_free(p_cg); + cfd_free(p_pcg); + cfd_free(p_temp); + cfd_free(rhs); +} + +/* ============================================================================ + * TEST: NON-2^K+1 DIMENSIONS ARE REJECTED AT INIT + * ============================================================================ */ + +void test_mg_pcg_rejects_invalid_dims(void) { + printf("\n Testing MG-PCG rejects non-2^k+1 grid dimensions...\n"); + + size_t n = 30; /* not 2^k+1 */ + double h = (DOMAIN_MAX - DOMAIN_MIN) / (n - 1); + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_CG, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL_MESSAGE(solver, "Could not create CG solver"); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.preconditioner = POISSON_PRECOND_MULTIGRID; + + cfd_status_t status = poisson_solver_init(solver, n, n, 1, h, h, 0.0, ¶ms); + TEST_ASSERT_EQUAL_INT_MESSAGE(CFD_ERROR_INVALID, status, + "CG init must propagate the inner MG rejection of non-2^k+1 dims"); + + poisson_solver_destroy(solver); +} + +/* ============================================================================ + * TEST: NON-SCALAR BACKENDS REJECT THE MG PRECONDITIONER + * ============================================================================ */ + +static void assert_backend_rejects_mg_precond(poisson_solver_method_t method, + poisson_solver_backend_t backend, + const char* label) { + poisson_solver_t* solver = poisson_solver_create(method, backend); + if (!solver) { + printf(" %s: not available (skipping)\n", label); + return; + } + + poisson_solver_params_t params = poisson_solver_params_default(); + params.preconditioner = POISSON_PRECOND_MULTIGRID; + + cfd_status_t status = poisson_solver_init(solver, 33, 33, 1, + 1.0 / 32.0, 1.0 / 32.0, 0.0, ¶ms); + printf(" %s: init status %d\n", label, (int)status); + TEST_ASSERT_EQUAL_INT_MESSAGE(CFD_ERROR_UNSUPPORTED, status, + "Backends without an MG preconditioner must reject it explicitly"); + + poisson_solver_destroy(solver); +} + +void test_mg_precond_unsupported_on_other_backends(void) { + printf("\n Testing MG preconditioner is rejected on unsupported backends...\n"); + + assert_backend_rejects_mg_precond(POISSON_METHOD_CG, POISSON_BACKEND_OMP, + "CG OMP"); + assert_backend_rejects_mg_precond(POISSON_METHOD_CG, POISSON_BACKEND_SIMD, + "CG SIMD"); + assert_backend_rejects_mg_precond(POISSON_METHOD_GMRES, POISSON_BACKEND_SCALAR, + "GMRES scalar"); + assert_backend_rejects_mg_precond(POISSON_METHOD_GMRES, POISSON_BACKEND_SIMD, + "GMRES SIMD"); +} + +/* ============================================================================ + * TEST: CONVENIENCE PRESET (POISSON_SOLVER_PCG_MG_SCALAR) + * ============================================================================ */ + +void test_pcg_mg_preset(void) { + printf("\n Testing poisson_solve_3d PCG_MG preset...\n"); + + /* Conforming grid: converges, returns iteration count */ + { + size_t n = 33; + double h = (DOMAIN_MAX - DOMAIN_MIN) / (n - 1); + double* p = create_field_3d(n, n, 1); + double* p_temp = create_field_3d(n, n, 1); + double* rhs = create_field_3d(n, n, 1); + TEST_ASSERT_NOT_NULL_MESSAGE(p, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_temp, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(rhs, "Memory allocation failed"); + + init_sinusoidal_rhs(rhs, n, n, h, h); + + int iters = poisson_solve_3d(p, p_temp, rhs, n, n, 1, h, h, 0.0, + POISSON_SOLVER_PCG_MG_SCALAR); + printf(" 33x33 preset solve: %d iterations\n", iters); + TEST_ASSERT_TRUE_MESSAGE(iters > 0, + "Preset must converge on a 2^k+1 grid and report iterations"); + + cfd_free(p); + cfd_free(p_temp); + cfd_free(rhs); + } + + /* Non-conforming grid: fails loudly with -1 (legacy contract) */ + { + size_t n = 30; + double h = (DOMAIN_MAX - DOMAIN_MIN) / (n - 1); + double* p = create_field_3d(n, n, 1); + double* p_temp = create_field_3d(n, n, 1); + double* rhs = create_field_3d(n, n, 1); + TEST_ASSERT_NOT_NULL_MESSAGE(p, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(p_temp, "Memory allocation failed"); + TEST_ASSERT_NOT_NULL_MESSAGE(rhs, "Memory allocation failed"); + + int iters = poisson_solve_3d(p, p_temp, rhs, n, n, 1, h, h, 0.0, + POISSON_SOLVER_PCG_MG_SCALAR); + printf(" 30x30 preset solve: returned %d\n", iters); + TEST_ASSERT_EQUAL_INT_MESSAGE(-1, iters, + "Preset must return -1 on non-2^k+1 grids"); + + cfd_free(p); + cfd_free(p_temp); + cfd_free(rhs); + } +} + +/* ============================================================================ + * MAIN + * ============================================================================ */ + +void setUp(void) {} +void tearDown(void) {} + +int main(void) { + UNITY_BEGIN(); + + printf("\n========================================\n"); + printf("MG-PCG Convergence Tests\n"); + printf("========================================\n"); + + RUN_TEST(test_mg_pcg_converges_correctly); + RUN_TEST(test_mg_pcg_iteration_reduction); + RUN_TEST(test_mg_pcg_3d); + RUN_TEST(test_mg_pcg_rejects_invalid_dims); + RUN_TEST(test_mg_precond_unsupported_on_other_backends); + RUN_TEST(test_pcg_mg_preset); + + printf("\n========================================\n"); + return UNITY_END(); +} From 96e0a29d8a707fda81423a6dc87176a587f46c91 Mon Sep 17 00:00:00 2001 From: shaia Date: Fri, 17 Jul 2026 17:50:30 +0300 Subject: [PATCH 2/9] Add pressure solver selection to the projection method New ns_solver_params_t.pressure_solver field selects the scalar projection's pressure Poisson solve: default CG (0, backward compatible), standalone multigrid V-cycles, or MG-preconditioned CG. The standalone mode subtracts the RHS interior mean first: the projection RHS is only near-compatible and the true Neumann system stalls on the incompatible component, while the CG-based presets are insensitive to it. Non-2^k+1 grids fail projection init with CFD_ERROR_UNSUPPORTED (probe-init of the MG hierarchy), and projection_optimized/omp/gpu reject any non-default selection at init: multigrid has no SIMD/OMP/GPU backend and cross-backend fallbacks are forbidden. --- CHANGELOG.md | 19 + CMakeLists.txt | 3 + README.md | 2 +- ROADMAP.md | 15 +- docs/reference/api-reference.md | 22 ++ docs/reference/solvers.md | 54 ++- .../cfd/solvers/navier_stokes_solver.h | 23 ++ lib/src/api/solver_registry.c | 49 ++- .../avx2/solver_projection_avx2.c | 7 +- .../navier_stokes/cpu/solver_explicit_euler.c | 3 +- .../navier_stokes/cpu/solver_projection.c | 30 +- .../cpu/test_projection_pressure_solver.c | 330 ++++++++++++++++++ 12 files changed, 543 insertions(+), 14 deletions(-) create mode 100644 tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c diff --git a/CHANGELOG.md b/CHANGELOG.md index 74eca05c..ad0cd3bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,25 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- **Multigrid wired into the projection method** — new + `ns_solver_params_t.pressure_solver` field (`ns_pressure_solver_t`; 0 = + existing CG behavior) selects the pressure Poisson solve of the scalar + `projection` solver: `NS_PRESSURE_SOLVER_MULTIGRID` (multigrid V-cycles, + RHS interior mean subtracted for Neumann compatibility) or + `NS_PRESSURE_SOLVER_PCG_MG` (MG-preconditioned CG). Non-2^k+1 grids and + non-scalar projection backends reject the selection with + `CFD_ERROR_UNSUPPORTED` at init + (`lib/src/solvers/navier_stokes/cpu/solver_projection.c`, + `lib/src/api/solver_registry.c`, + `tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c`). +- **Multigrid-preconditioned CG** — `POISSON_PRECOND_MULTIGRID` runs one + symmetric V(2,2) weighted-Jacobi multigrid cycle in Dirichlet mode per + preconditioner apply (scalar CG only; other CG/GMRES backends reject it + with `CFD_ERROR_UNSUPPORTED`). Grid-size-independent convergence: 5 CG + iterations at 33²–129² (tol 1e-8) vs 50–170 unpreconditioned. Exposed as + the `POISSON_SOLVER_PCG_MG_SCALAR` convenience preset + (`lib/src/solvers/linear/cpu/linear_solver_cg.c`, + `tests/math/test_mg_pcg_convergence.c`). - **Geometric multigrid Poisson solver** (scalar backend) — V/W/F(FMG) cycles with Red-Black Gauss-Seidel or weighted-Jacobi smoothers, full-weighting restriction and bilinear/trilinear prolongation, 2D/3D. Two BC modes: Neumann zero-gradient (default, diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d161b7f..a6365a15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -253,6 +253,7 @@ if(BUILD_TESTS) add_executable(test_solver_rk2 tests/solvers/navier_stokes/cpu/test_solver_rk2.c) add_executable(test_solver_rk4 tests/solvers/navier_stokes/cpu/test_solver_rk4.c) add_executable(test_solver_projection tests/solvers/navier_stokes/cpu/test_solver_projection.c) + add_executable(test_projection_pressure_solver tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c) add_executable(test_ns_solver_3d tests/solvers/navier_stokes/cpu/test_ns_solver_3d.c) # Navier-Stokes solver tests (AVX2) @@ -408,6 +409,7 @@ if(BUILD_TESTS) target_link_libraries(test_solver_rk2 PRIVATE CFD::Library unity) target_link_libraries(test_solver_rk4 PRIVATE CFD::Library unity) target_link_libraries(test_solver_projection PRIVATE CFD::Library unity) + target_link_libraries(test_projection_pressure_solver PRIVATE CFD::Library unity) target_link_libraries(test_ns_solver_3d PRIVATE CFD::Library unity) target_link_libraries(test_solver_explicit_euler_avx2 PRIVATE CFD::Library unity) target_link_libraries(test_solver_projection_avx2 PRIVATE CFD::Library unity) @@ -523,6 +525,7 @@ if(BUILD_TESTS) add_test(NAME SolverRK2Test COMMAND test_solver_rk2) add_test(NAME SolverRK4Test COMMAND test_solver_rk4) add_test(NAME SolverProjectionTest COMMAND test_solver_projection) + add_test(NAME ProjectionPressureSolverTest COMMAND test_projection_pressure_solver) add_test(NAME NSSolver3DTest COMMAND test_ns_solver_3d) add_test(NAME SolverExplicitEulerAvx2Test COMMAND test_solver_explicit_euler_avx2) add_test(NAME SolverProjectionAvx2Test COMMAND test_solver_projection_avx2) diff --git a/README.md b/README.md index 09bdfa69..3c7d8063 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ A production-grade computational fluid dynamics (CFD) library in C for solving 2 - 🚀 **Multiple Backends**: CPU (scalar), SIMD (AVX2/NEON), OpenMP, CUDA - 🔧 **Pluggable Solvers**: Explicit Euler, RK2 (Heun), Projection Method (Chorin's algorithm) -- 📊 **Linear Solvers**: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB +- 📊 **Linear Solvers**: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB, GMRES(m), geometric multigrid (V/W/F cycles, also as CG preconditioner and projection pressure solver via `params.pressure_solver`) - 🌡️ **Heat Transfer**: Energy equation (advection–diffusion) + Boussinesq buoyancy + thermal BCs - 🌀 **Turbulence (RANS)**: Standard k-ε and Spalart-Allmaras with log-law wall functions; validated vs turbulent channel flow at Re_τ = 395 - 🎯 **Validated**: Ghia lid-driven cavity, Taylor-Green vortex, Poiseuille flow, natural convection, turbulent channel flow diff --git a/ROADMAP.md b/ROADMAP.md index 587d1a50..88a57a38 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -76,7 +76,9 @@ Genuine constraints to be aware of (not backlog items): Windows/macOS handle automatically. Future option: weak symbols or a plugin architecture. - **Multigrid grid-dimension constraint** — the geometric multigrid solver requires 2^k+1 points per active dimension (e.g. 33, 65, 129) for exact coarsening; other sizes - return `CFD_ERROR_INVALID`. Use CG for arbitrary grid sizes. + return `CFD_ERROR_INVALID`. Use CG for arbitrary grid sizes. On conforming grids, + MG-preconditioned CG (`POISSON_PRECOND_MULTIGRID`) gives grid-size-independent + iteration counts on uniform grids where Jacobi PCG gives no benefit. > The SOR optimal-ω and Jacobi spectral-radius formulas (ρ = cos(πh), > ω = 2/(1 + sin(πh))) assume Dirichlet BCs; with Neumann BCs optimal ω is typically lower @@ -130,8 +132,10 @@ non-symmetric operators arrive, e.g. implicit advection-diffusion in §1.5). - [x] Geometric multigrid — scalar backend; V/W/F(FMG) cycles, Red-Black GS or weighted Jacobi smoothers, full-weighting restriction (Neumann-folded at boundaries) + bilinear/trilinear prolongation, Neumann (default) and Dirichlet BC modes, 2D/3D. - Grid dims must be 2^k+1. SIMD/OMP/GPU variants and MG-preconditioned CG deferred - (see `.claude/specs/multigrid-projection-integration.md` for projection wiring). + Grid dims must be 2^k+1. Wired into the projection method + (`ns_solver_params_t.pressure_solver`, scalar backend) and available as a CG + preconditioner (`POISSON_PRECOND_MULTIGRID`, scalar CG). SIMD/OMP/GPU variants + deferred (see `.claude/specs/multigrid-projection-integration.md`). - [ ] Algebraic multigrid (AMG) — solver and preconditioner (for CG/GMRES/BiCGSTAB) - [x] GPU plain SOR (Block SOR: per-thread tile sweep, red-black tile coloring, in-place; closes the matrix) @@ -209,7 +213,10 @@ SIMD Poisson integration is done; current ~1.3–1.5× speedup is Amdahl-limited (parallelizable fraction ~80%). Remaining optimization work: - [ ] Increase `POISSON_MAX_ITER` or implement adaptive tolerance -- [ ] Optional multigrid preconditioner for faster convergence (see §1.2 multigrid) +- [x] Optional multigrid preconditioner for faster convergence (see §1.2 multigrid) — + done for scalar CG (`POISSON_PRECOND_MULTIGRID`, symmetric V(2,2) Jacobi cycle) and + selectable in the scalar projection via `pressure_solver = NS_PRESSURE_SOLVER_PCG_MG`; + SIMD/OMP variants follow the multigrid backend work in §1.2 - [ ] Red-Black omega parameter tuning - [ ] Profile to identify remaining bottlenecks - [ ] OpenMP+SIMD hybrid projection (OMP across rows, SIMD within rows); benchmark vs pure OMP diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md index b25484ed..17a2a4e4 100644 --- a/docs/reference/api-reference.md +++ b/docs/reference/api-reference.md @@ -263,11 +263,26 @@ typedef struct { // Turbulence turbulence_model_t turb_model; // TURB_MODEL_NONE (default), K_EPSILON, SPALART_ALLMARAS ns_turbulence_bc_config_t turb_bc; // Per-face turbulence BC types + // Pressure Poisson solver (projection solvers) + ns_pressure_solver_t pressure_solver; // NS_PRESSURE_SOLVER_DEFAULT (0) = backend's CG } ns_solver_params_t; ns_solver_params_t ns_solver_params_default(void); ``` +`ns_pressure_solver_t` selects the pressure Poisson solve of the projection +method (scalar `projection` solver only; the multigrid modes require 2^k+1 +grid points per active dimension, and other projection backends reject +non-default values with `CFD_ERROR_UNSUPPORTED` at init): + +```c +typedef enum { + NS_PRESSURE_SOLVER_DEFAULT = 0, // Backend's default CG pressure solve + NS_PRESSURE_SOLVER_MULTIGRID = 1, // Geometric multigrid V-cycle solver + NS_PRESSURE_SOLVER_PCG_MG = 2, // CG preconditioned by one MG V-cycle +} ns_pressure_solver_t; +``` + `turbulence_model_t` is defined in `cfd/solvers/navier_stokes_solver.h`: ```c @@ -492,6 +507,13 @@ poisson_solver_params_t poisson_solver_params_default(void); > `poisson_solver_init` returns `CFD_ERROR_INVALID` otherwise. In the default > `MG_BC_NEUMANN` mode the solution is defined up to an additive constant. +`poisson_precond_type_t` values: `POISSON_PRECOND_NONE` (0, default), +`POISSON_PRECOND_JACOBI` (1), `POISSON_PRECOND_MULTIGRID` (2 — one MG V-cycle +per apply; scalar CG only, 2^k+1 dims; other backends return +`CFD_ERROR_UNSUPPORTED`). The convenience API exposes the MG-preconditioned CG +as the `POISSON_SOLVER_PCG_MG_SCALAR` preset for +`poisson_solve()`/`poisson_solve_3d()`. + ### Poisson Statistics ```c diff --git a/docs/reference/solvers.md b/docs/reference/solvers.md index 037a944b..6e4978a9 100644 --- a/docs/reference/solvers.md +++ b/docs/reference/solvers.md @@ -78,6 +78,30 @@ Chorin's projection method - properly enforces incompressibility constraint. | `projection_omp` | OpenMP | Multi-threaded | | `projection_jacobi_gpu` | GPU | CUDA-accelerated (Jacobi iteration) | +**Pressure solver selection** (`ns_solver_params_t.pressure_solver`): + +Each backend pairs with a CG Poisson preset by default. On the scalar +`projection` solver the pressure solve can be switched to geometric multigrid: + +| `pressure_solver` value | Pressure Poisson solve | +|-------------------------|------------------------| +| `NS_PRESSURE_SOLVER_DEFAULT` (0) | Backend's CG preset (existing behavior) | +| `NS_PRESSURE_SOLVER_MULTIGRID` | Multigrid V-cycles (`POISSON_SOLVER_MG_SCALAR`) | +| `NS_PRESSURE_SOLVER_PCG_MG` | CG + MG V-cycle preconditioner (`POISSON_SOLVER_PCG_MG_SCALAR`) | + +The MG modes require 2^k+1 grid points per active dimension (33, 65, 129, ...); +`solver_init` returns `CFD_ERROR_UNSUPPORTED` otherwise. `projection_optimized`, +`projection_omp`, and `projection_gpu` reject any non-default value with +`CFD_ERROR_UNSUPPORTED` at init — multigrid has no SIMD/OMP/GPU backend yet and +the library never falls back across backends silently. + +```c +ns_solver_params_t params = ns_solver_params_default(); +params.pressure_solver = NS_PRESSURE_SOLVER_MULTIGRID; // 2^k+1 grids only +ns_solver_t* slv = cfd_solver_create(registry, NS_SOLVER_TYPE_PROJECTION); +cfd_status_t status = solver_init(slv, grid, ¶ms); // UNSUPPORTED on 128x128 +``` + ## Linear Solvers (Poisson Equation) The projection method requires solving the pressure Poisson equation: @@ -196,6 +220,16 @@ poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_CG, - Simple, cheap per iteration - Variable coefficients: can reduce iterations - Constant coefficients (uniform grid): no benefit +- **Multigrid** (`POISSON_PRECOND_MULTIGRID`): one geometric-multigrid V-cycle + per apply + - Grid-size-independent outer iteration count (measured 5 CG iterations at + 33²–129², tol 1e-8, vs 50–170 for plain CG) + - Scalar CG backend only; other backends reject it with + `CFD_ERROR_UNSUPPORTED` + - Requires 2^k+1 grid points per active dimension (inherited from the + multigrid hierarchy); init fails with `CFD_ERROR_INVALID` otherwise + - The inner cycle is a symmetric V(2,2) with weighted-Jacobi smoothing in + Dirichlet mode, so the preconditioner is SPD as CG requires - **SSOR**: M = (D + L)D⁻¹(D + U) (future) - **ILU**: Incomplete LU factorization (future) @@ -209,6 +243,19 @@ poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_CG, poisson_solver_init(solver, nx, ny, dx, dy, ¶ms); // Pass params with preconditioner ``` +**Multigrid-preconditioned CG (scalar backend, 2^k+1 dims):** + +```c +poisson_solver_params_t params = poisson_solver_params_default(); +params.preconditioner = POISSON_PRECOND_MULTIGRID; + +poisson_solver_t* solver = poisson_solver_create(POISSON_METHOD_CG, + POISSON_BACKEND_SCALAR); +poisson_solver_init(solver, 65, 65, 1, dx, dy, 0.0, ¶ms); +``` +The convenience API exposes the same configuration as the +`POISSON_SOLVER_PCG_MG_SCALAR` preset for `poisson_solve()`/`poisson_solve_3d()`. + #### 6. BiCGSTAB **Algorithm:** Bi-Conjugate Gradient Stabilized for non-symmetric systems. @@ -282,8 +329,10 @@ Smoothers: Red-Black Gauss-Seidel (default) or weighted Jacobi (ω=2/3). - Parameters: `mg_cycle`, `mg_smoother`, `mg_bc`, `mg_pre_smooth`/`mg_post_smooth` (default 2/2), `mg_coarse_max_iter` (default 50), `mg_max_levels` (0 = auto) -**Backends:** scalar only. SIMD/OMP/GPU variants and use as a CG preconditioner -are planned follow-ups. +**Backends:** scalar only. Also available as a CG preconditioner +(`POISSON_PRECOND_MULTIGRID`, scalar CG only — see §5) and as the projection +method's pressure solver (`ns_solver_params_t.pressure_solver`). SIMD/OMP/GPU +variants are planned follow-ups. **Usage:** ```c @@ -307,6 +356,7 @@ poisson_solver_init(solver, 65, 65, 1, dx, dy, 0.0, ¶ms); // dims 2^k+1 | Red-Black SOR | ~2000 | 8 | Parallelizable | | CG | ~80 | 5 | Best for large grids | | PCG (Jacobi) | ~80 | 5.5 | No benefit on uniform grid | +| PCG (Multigrid) | ~5 | — | Grid-size-independent; needs 2^k+1 dims | | BiCGSTAB | ~40 | 4 | Fastest convergence | | Multigrid V(2,2) | ~8 cycles | — | O(N), grid-size-independent; needs 2^k+1 dims | diff --git a/lib/include/cfd/solvers/navier_stokes_solver.h b/lib/include/cfd/solvers/navier_stokes_solver.h index 6e12b869..463a0bcc 100644 --- a/lib/include/cfd/solvers/navier_stokes_solver.h +++ b/lib/include/cfd/solvers/navier_stokes_solver.h @@ -160,6 +160,25 @@ typedef struct { bc_dirichlet_values_t nu_tilde_values; /**< Fixed nu_tilde per Dirichlet face */ } ns_turbulence_bc_config_t; +/** + * Pressure Poisson solver selection for projection-method solvers. + * + * NS_PRESSURE_SOLVER_DEFAULT (0) keeps each backend's existing CG pressure + * solve, so zero-initialization is fully backward compatible. + * + * The multigrid modes are implemented on the scalar "projection" solver only + * and require 2^k+1 grid points per active dimension (e.g. 33, 65, 129). + * projection_optimized, projection_omp, and projection_gpu reject any + * non-default value with CFD_ERROR_UNSUPPORTED at init (no silent + * cross-backend fallbacks); non-conforming grid dimensions are likewise + * rejected at init with CFD_ERROR_UNSUPPORTED. + */ +typedef enum { + NS_PRESSURE_SOLVER_DEFAULT = 0, /**< Backend's default CG pressure solve */ + NS_PRESSURE_SOLVER_MULTIGRID = 1, /**< Geometric multigrid V-cycle solver */ + NS_PRESSURE_SOLVER_PCG_MG = 2, /**< CG preconditioned by one MG V-cycle */ +} ns_pressure_solver_t; + /** * Navier-Stokes solver parameters */ @@ -207,6 +226,10 @@ typedef struct { * a turbulence model is enabled. */ turbulence_model_t turb_model; /**< Turbulence model selection */ ns_turbulence_bc_config_t turb_bc; /**< Per-face turbulence BCs (zero-init = all PERIODIC) */ + + /* Pressure Poisson solver for projection solvers (scalar "projection" + * backend only; 0 = existing CG behavior, backward compatible). */ + ns_pressure_solver_t pressure_solver; /**< Pressure solver selection */ } ns_solver_params_t; diff --git a/lib/src/api/solver_registry.c b/lib/src/api/solver_registry.c index e0dff077..2d6d437d 100644 --- a/lib/src/api/solver_registry.c +++ b/lib/src/api/solver_registry.c @@ -919,8 +919,32 @@ typedef struct { } projection_context; static cfd_status_t projection_init(ns_solver_t* solver, const grid* grid, const ns_solver_params_t* params) { - (void)grid; - (void)params; + if (!solver || !grid) { + return CFD_ERROR_INVALID; + } + if (params && params->pressure_solver != NS_PRESSURE_SOLVER_DEFAULT) { + if (params->pressure_solver != NS_PRESSURE_SOLVER_MULTIGRID && + params->pressure_solver != NS_PRESSURE_SOLVER_PCG_MG) { + return CFD_ERROR_INVALID; + } + /* Both MG modes require a multigrid hierarchy on this exact grid: + * probe-init and reject non-2^k+1 dimensions up front. */ + poisson_solver_t* probe = poisson_solver_create( + POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); + if (!probe) { + return CFD_ERROR_UNSUPPORTED; + } + cfd_status_t probe_status = poisson_solver_init( + probe, grid->nx, grid->ny, grid->nz, + grid->dx[0], grid->dy[0], + (grid->nz > 1 && grid->dz) ? grid->dz[0] : 0.0, NULL); + poisson_solver_destroy(probe); + if (probe_status != CFD_SUCCESS) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, + "Multigrid pressure solver requires 2^k+1 grid points per active dimension"); + return CFD_ERROR_UNSUPPORTED; + } + } projection_context* ctx = (projection_context*)cfd_malloc(sizeof(projection_context)); if (!ctx) { return CFD_ERROR; @@ -1142,6 +1166,18 @@ static ns_solver_t* create_explicit_euler_gpu_solver(void) { * Built-in solver: GPU-Accelerated Projection Method */ +static cfd_status_t gpu_projection_init(ns_solver_t* solver, const grid* grid, + const ns_solver_params_t* params) { + (void)solver; + (void)grid; + if (params && params->pressure_solver != NS_PRESSURE_SOLVER_DEFAULT) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, + "Multigrid pressure solver is only supported by the scalar projection solver"); + return CFD_ERROR_UNSUPPORTED; + } + return CFD_SUCCESS; +} + static cfd_status_t gpu_projection_step(ns_solver_t* solver, flow_field* field, const grid* grid, const ns_solver_params_t* params, ns_solver_stats_t* stats) { (void)solver; @@ -1196,7 +1232,7 @@ static ns_solver_t* create_projection_gpu_solver(void) { s->capabilities = NS_SOLVER_CAP_INCOMPRESSIBLE | NS_SOLVER_CAP_TRANSIENT | NS_SOLVER_CAP_GPU; s->backend = NS_SOLVER_BACKEND_CUDA; - s->init = NULL; + s->init = gpu_projection_init; // Rejects non-default pressure_solver s->destroy = NULL; s->step = gpu_projection_step; s->solve = gpu_projection_solve; @@ -1535,11 +1571,16 @@ static ns_solver_t* create_rk4_omp_solver(void) { static cfd_status_t projection_omp_init(ns_solver_t* solver, const grid* grid, const ns_solver_params_t* params) { - (void)params; if (!solver || !grid) { return CFD_ERROR_INVALID; } + if (params && params->pressure_solver != NS_PRESSURE_SOLVER_DEFAULT) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, + "Multigrid pressure solver is only supported by the scalar projection solver"); + return CFD_ERROR_UNSUPPORTED; + } + /* OMP projection requires OMP CG Poisson solver. * Never fall back to scalar CG — it would serialize the Poisson solve. */ poisson_solver_t* test_solver = poisson_solver_create( diff --git a/lib/src/solvers/navier_stokes/avx2/solver_projection_avx2.c b/lib/src/solvers/navier_stokes/avx2/solver_projection_avx2.c index e9412ee0..7dc0ac9e 100644 --- a/lib/src/solvers/navier_stokes/avx2/solver_projection_avx2.c +++ b/lib/src/solvers/navier_stokes/avx2/solver_projection_avx2.c @@ -80,7 +80,6 @@ cfd_status_t projection_simd_step(struct NSSolver* solver, flow_field* field, co cfd_status_t projection_simd_init(struct NSSolver* solver, const grid* grid, const ns_solver_params_t* params) { - (void)params; if (!solver || !grid) { return CFD_ERROR_INVALID; } @@ -88,6 +87,12 @@ cfd_status_t projection_simd_init(struct NSSolver* solver, const grid* grid, return CFD_ERROR_INVALID; } + if (params && params->pressure_solver != NS_PRESSURE_SOLVER_DEFAULT) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, + "Multigrid pressure solver is only supported by the scalar projection solver"); + return CFD_ERROR_UNSUPPORTED; + } + /* Verify SIMD CG Poisson solver is available before allocating resources */ poisson_solver_t* test_solver = poisson_solver_create( POISSON_METHOD_CG, POISSON_BACKEND_SIMD); diff --git a/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c b/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c index 7aa9320e..d8123fed 100644 --- a/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c +++ b/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c @@ -77,7 +77,8 @@ ns_solver_params_t ns_solver_params_default(void) { .heat_source_context = NULL, .thermal_bc = {0}, .turb_model = TURB_MODEL_NONE, - .turb_bc = {0}}; + .turb_bc = {0}, + .pressure_solver = NS_PRESSURE_SOLVER_DEFAULT}; return params; } flow_field* flow_field_create(size_t nx, size_t ny, size_t nz) { diff --git a/lib/src/solvers/navier_stokes/cpu/solver_projection.c b/lib/src/solvers/navier_stokes/cpu/solver_projection.c index b479dda5..04334c56 100644 --- a/lib/src/solvers/navier_stokes/cpu/solver_projection.c +++ b/lib/src/solvers/navier_stokes/cpu/solver_projection.c @@ -27,6 +27,7 @@ #include "cfd/solvers/turbulence_solver.h" #include "../../energy/energy_solver_internal.h" +#include "../../linear/multigrid_internal.h" #include "../../turbulence/turbulence_solver_internal.h" #include "../boundary_copy_utils.h" @@ -78,6 +79,24 @@ cfd_status_t solve_projection_method(flow_field* field, const grid* grid, double dt = params->dt; double nu = params->mu; + /* Map the pressure-solver selection to a Poisson preset. Grid-dimension + * compatibility for the MG modes is validated at solver init; a failed + * solve here still degrades loudly via poisson_iters < 0. */ + poisson_solver_type pressure_preset; + switch (params->pressure_solver) { + case NS_PRESSURE_SOLVER_DEFAULT: + pressure_preset = POISSON_SOLVER_CG_SCALAR; + break; + case NS_PRESSURE_SOLVER_MULTIGRID: + pressure_preset = POISSON_SOLVER_MG_SCALAR; + break; + case NS_PRESSURE_SOLVER_PCG_MG: + pressure_preset = POISSON_SOLVER_PCG_MG_SCALAR; + break; + default: + return CFD_ERROR_INVALID; + } + /* Branch-free 3D constants */ size_t stride_z = (nz > 1) ? plane : 0; size_t k_start = (nz > 1) ? 1 : 0; @@ -250,9 +269,18 @@ cfd_status_t solve_projection_method(flow_field* field, const grid* grid, } } + /* Neumann compatibility projection: standalone multigrid solves the + * true singular Neumann system, so the RHS must have zero interior + * mean or the residual stalls at the incompatible component. The + * CG-based presets are insensitive to it (their interior-only Krylov + * updates act as a nonsingular operator) and keep today's behavior. */ + if (pressure_preset == POISSON_SOLVER_MG_SCALAR) { + mg_subtract_interior_mean(rhs, nx, ny, nz); + } + /* Solve Poisson equation using library solver */ int poisson_iters = poisson_solve_3d(p_new, p_temp, rhs, nx, ny, nz, dx, dy, dz, - POISSON_SOLVER_CG_SCALAR); + pressure_preset); if (poisson_iters < 0) { cfd_free(u_star); cfd_free(v_star); cfd_free(w_star); diff --git a/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c b/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c new file mode 100644 index 00000000..aeeb3a8b --- /dev/null +++ b/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c @@ -0,0 +1,330 @@ +/** + * Pressure-solver selection tests for the projection method + * + * Verifies the ns_solver_params_t.pressure_solver field: + * - NS_PRESSURE_SOLVER_MULTIGRID and NS_PRESSURE_SOLVER_PCG_MG produce + * divergence-free fields on the scalar "projection" solver + * - The MG modes agree with the default CG pressure solve within solver + * tolerance + * - Non-2^k+1 grids are rejected at init with CFD_ERROR_UNSUPPORTED + * - projection_optimized / projection_omp / projection_gpu reject any + * non-default selection with CFD_ERROR_UNSUPPORTED (no silent fallbacks) + * - Zero-initialized params keep the existing CG behavior bit-for-bit + */ + +#include "../test_solver_helpers.h" +#include "cfd/core/cfd_init.h" +#include "cfd/core/cfd_status.h" +#include "cfd/core/grid.h" +#include "cfd/core/indexing.h" +#include "cfd/core/memory.h" +#include "cfd/solvers/navier_stokes_solver.h" +#include "unity.h" + +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +void setUp(void) { + cfd_init(); +} + +void tearDown(void) { + cfd_finalize(); +} + +#define GRID_N 33 /* 2^5+1: multigrid-conforming */ +#define GRID_N_BAD 30 /* not 2^k+1 */ +#define NUM_STEPS 5 +#define TEST_DT 1e-3 + +static ns_solver_params_t make_params(ns_pressure_solver_t pressure_solver) { + ns_solver_params_t params = ns_solver_params_default(); + params.dt = TEST_DT; + params.mu = 0.01; + params.max_iter = 1; + params.pressure_solver = pressure_solver; + return params; +} + +/** + * Run the scalar "projection" solver for NUM_STEPS on a Taylor-Green field. + * Returns the init status; on CFD_SUCCESS the field holds the final state. + */ +static cfd_status_t run_projection_taylor_green(ns_pressure_solver_t pressure_solver, + flow_field* field, const grid* g) { + ns_solver_params_t params = make_params(pressure_solver); + + ns_solver_registry_t* registry = cfd_registry_create(); + TEST_ASSERT_NOT_NULL_MESSAGE(registry, "Failed to create registry"); + cfd_registry_register_defaults(registry); + + ns_solver_t* slv = cfd_solver_create(registry, NS_SOLVER_TYPE_PROJECTION); + TEST_ASSERT_NOT_NULL_MESSAGE(slv, "projection solver not available"); + + cfd_status_t init_status = solver_init(slv, g, ¶ms); + if (init_status != CFD_SUCCESS) { + solver_destroy(slv); + cfd_registry_destroy(registry); + return init_status; + } + + test_init_taylor_green(field, g); + + ns_solver_stats_t stats = ns_solver_stats_default(); + for (int step = 0; step < NUM_STEPS; step++) { + cfd_status_t step_status = solver_step(slv, field, g, ¶ms, &stats); + TEST_ASSERT_EQUAL_INT_MESSAGE(CFD_SUCCESS, step_status, + "projection step failed"); + } + + solver_destroy(slv); + cfd_registry_destroy(registry); + return CFD_SUCCESS; +} + +//============================================================================= +// TEST: MG PRESSURE SOLVES PRODUCE DIVERGENCE-FREE FIELDS +//============================================================================= + +static void run_divergence_free_check(ns_pressure_solver_t pressure_solver, + const char* label) { + ns_solver_params_t params = make_params(pressure_solver); + + /* Divergence tolerance is relaxed (matching test_solver_projection.c): + * the projection solver's simple iterative method may not fully converge + * in a few steps. The meaningful assertion is that each step reduces + * the divergence. */ + test_result result = test_run_divergence_free( + NS_SOLVER_TYPE_PROJECTION, GRID_N, GRID_N, ¶ms, 10, 1.0); + + if (result.solver_unavailable) { + printf(" %s: solver unavailable, skipping\n", label); + return; + } + + printf(" %s: %s\n", label, result.message); + TEST_ASSERT_TRUE_MESSAGE(result.passed, result.message); + TEST_ASSERT_TRUE_MESSAGE(result.final_divergence < result.initial_divergence, + "Projection with an MG pressure solve must reduce divergence"); +} + +void test_projection_mg_divergence_free(void) { + printf("\n Testing projection + multigrid pressure solve...\n"); + run_divergence_free_check(NS_PRESSURE_SOLVER_MULTIGRID, "MULTIGRID"); +} + +void test_projection_pcg_mg_divergence_free(void) { + printf("\n Testing projection + MG-preconditioned CG pressure solve...\n"); + run_divergence_free_check(NS_PRESSURE_SOLVER_PCG_MG, "PCG_MG"); +} + +//============================================================================= +// TEST: MG PRESSURE SOLVES AGREE WITH THE DEFAULT CG SOLVE +//============================================================================= + +void test_projection_mg_matches_cg(void) { + printf("\n Testing MG pressure solves agree with default CG...\n"); + + grid* g = grid_create(GRID_N, GRID_N, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + flow_field* field_cg = flow_field_create(GRID_N, GRID_N, 1); + flow_field* field_mg = flow_field_create(GRID_N, GRID_N, 1); + flow_field* field_pcg = flow_field_create(GRID_N, GRID_N, 1); + TEST_ASSERT_NOT_NULL(field_cg); + TEST_ASSERT_NOT_NULL(field_mg); + TEST_ASSERT_NOT_NULL(field_pcg); + + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, + run_projection_taylor_green(NS_PRESSURE_SOLVER_DEFAULT, field_cg, g)); + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, + run_projection_taylor_green(NS_PRESSURE_SOLVER_MULTIGRID, field_mg, g)); + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, + run_projection_taylor_green(NS_PRESSURE_SOLVER_PCG_MG, field_pcg, g)); + + size_t n = (size_t)GRID_N * GRID_N; + double u_norm = test_compute_l2_norm(field_cg->u, n); + double v_norm = test_compute_l2_norm(field_cg->v, n); + TEST_ASSERT_TRUE_MESSAGE(u_norm > 0.0 && v_norm > 0.0, + "CG reference run produced a zero field"); + + double u_err_mg = test_compute_l2_error(field_cg->u, field_mg->u, n) / u_norm; + double v_err_mg = test_compute_l2_error(field_cg->v, field_mg->v, n) / v_norm; + double u_err_pcg = test_compute_l2_error(field_cg->u, field_pcg->u, n) / u_norm; + double v_err_pcg = test_compute_l2_error(field_cg->v, field_pcg->v, n) / v_norm; + + printf(" MULTIGRID vs CG: rel L2 u=%.2e v=%.2e\n", u_err_mg, v_err_mg); + printf(" PCG_MG vs CG: rel L2 u=%.2e v=%.2e\n", u_err_pcg, v_err_pcg); + + /* Pressure is skipped: Neumann solutions differ by an additive constant. + * + * PCG_MG shares CG's effective operator (interior-only Krylov updates), + * so it matches CG almost exactly. Standalone multigrid solves the true + * Neumann system with a mean-compatible RHS — a slightly different + * boundary treatment of the pressure solve — so velocity agreement is + * O(1e-3) after a few steps (measured 3e-3 on 33x33), not bitwise. */ + TEST_ASSERT_TRUE_MESSAGE(u_err_mg < TOLERANCE_RELAXED && v_err_mg < TOLERANCE_RELAXED, + "Multigrid pressure solve should agree with CG within solver tolerance"); + TEST_ASSERT_TRUE_MESSAGE(u_err_pcg < TOLERANCE_MODERATE && v_err_pcg < TOLERANCE_MODERATE, + "MG-PCG pressure solve should agree with CG within solver tolerance"); + + flow_field_destroy(field_cg); + flow_field_destroy(field_mg); + flow_field_destroy(field_pcg); + grid_destroy(g); +} + +//============================================================================= +// TEST: NON-2^K+1 GRID IS REJECTED AT INIT +//============================================================================= + +void test_projection_mg_rejects_non_pow2_grid(void) { + printf("\n Testing MG pressure solver rejects non-2^k+1 grids at init...\n"); + + grid* g = grid_create(GRID_N_BAD, GRID_N_BAD, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + ns_solver_registry_t* registry = cfd_registry_create(); + TEST_ASSERT_NOT_NULL(registry); + cfd_registry_register_defaults(registry); + + ns_pressure_solver_t modes[] = { + NS_PRESSURE_SOLVER_MULTIGRID, + NS_PRESSURE_SOLVER_PCG_MG, + }; + for (size_t m = 0; m < sizeof(modes) / sizeof(modes[0]); m++) { + ns_solver_t* slv = cfd_solver_create(registry, NS_SOLVER_TYPE_PROJECTION); + TEST_ASSERT_NOT_NULL_MESSAGE(slv, "projection solver not available"); + + ns_solver_params_t params = make_params(modes[m]); + cfd_status_t init_status = solver_init(slv, g, ¶ms); + printf(" mode %d on %dx%d: init status %d\n", + (int)modes[m], GRID_N_BAD, GRID_N_BAD, (int)init_status); + TEST_ASSERT_EQUAL_INT_MESSAGE(CFD_ERROR_UNSUPPORTED, init_status, + "MG pressure solver on a non-2^k+1 grid must fail init with UNSUPPORTED"); + + solver_destroy(slv); + } + + cfd_registry_destroy(registry); + grid_destroy(g); +} + +//============================================================================= +// TEST: NON-SCALAR PROJECTION BACKENDS REJECT MG SELECTION +//============================================================================= + +void test_projection_backends_reject_mg(void) { + printf("\n Testing non-scalar projection backends reject MG selection...\n"); + + /* Conforming 33x33 grid: proves the rejection is an explicit policy, + * not a dimensional failure. */ + grid* g = grid_create(GRID_N, GRID_N, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + ns_solver_registry_t* registry = cfd_registry_create(); + TEST_ASSERT_NOT_NULL(registry); + cfd_registry_register_defaults(registry); + + const char* backends[] = { + NS_SOLVER_TYPE_PROJECTION_OPTIMIZED, + NS_SOLVER_TYPE_PROJECTION_OMP, + NS_SOLVER_TYPE_PROJECTION_GPU, + }; + ns_pressure_solver_t modes[] = { + NS_PRESSURE_SOLVER_MULTIGRID, + NS_PRESSURE_SOLVER_PCG_MG, + }; + + for (size_t b = 0; b < sizeof(backends) / sizeof(backends[0]); b++) { + for (size_t m = 0; m < sizeof(modes) / sizeof(modes[0]); m++) { + ns_solver_t* slv = cfd_solver_create(registry, backends[b]); + if (!slv) { + printf(" %s: not available (skipping)\n", backends[b]); + break; + } + + ns_solver_params_t params = make_params(modes[m]); + cfd_status_t init_status = solver_init(slv, g, ¶ms); + printf(" %s, mode %d: init status %d\n", + backends[b], (int)modes[m], (int)init_status); + TEST_ASSERT_EQUAL_INT_MESSAGE(CFD_ERROR_UNSUPPORTED, init_status, + "Non-scalar projection backends must reject MG pressure solvers"); + + solver_destroy(slv); + } + } + + cfd_registry_destroy(registry); + grid_destroy(g); +} + +//============================================================================= +// TEST: ZERO-INITIALIZED PARAMS KEEP THE EXISTING CG BEHAVIOR +//============================================================================= + +void test_projection_zero_init_backward_compat(void) { + printf("\n Testing zero-init params preserve default CG behavior...\n"); + + ns_solver_params_t defaults = ns_solver_params_default(); + TEST_ASSERT_EQUAL_INT_MESSAGE(NS_PRESSURE_SOLVER_DEFAULT, defaults.pressure_solver, + "ns_solver_params_default must leave pressure_solver at the zero default"); + + grid* g = grid_create(GRID_N, GRID_N, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + flow_field* field_a = flow_field_create(GRID_N, GRID_N, 1); + flow_field* field_b = flow_field_create(GRID_N, GRID_N, 1); + TEST_ASSERT_NOT_NULL(field_a); + TEST_ASSERT_NOT_NULL(field_b); + + /* Same run twice: once relying on the default, once with the explicit + * enum value. The scalar path is deterministic, so bitwise equality. */ + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, + run_projection_taylor_green(NS_PRESSURE_SOLVER_DEFAULT, field_a, g)); + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, + run_projection_taylor_green((ns_pressure_solver_t)0, field_b, g)); + + size_t bytes = (size_t)GRID_N * GRID_N * sizeof(double); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_a->u, field_b->u, bytes), + "u fields must be bitwise identical"); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_a->v, field_b->v, bytes), + "v fields must be bitwise identical"); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_a->p, field_b->p, bytes), + "p fields must be bitwise identical"); + + flow_field_destroy(field_a); + flow_field_destroy(field_b); + grid_destroy(g); +} + +//============================================================================= +// MAIN +//============================================================================= + +int main(void) { + UNITY_BEGIN(); + + printf("\n========================================\n"); + printf("Projection Pressure-Solver Selection Tests\n"); + printf("========================================\n"); + + RUN_TEST(test_projection_mg_divergence_free); + RUN_TEST(test_projection_pcg_mg_divergence_free); + RUN_TEST(test_projection_mg_matches_cg); + RUN_TEST(test_projection_mg_rejects_non_pow2_grid); + RUN_TEST(test_projection_backends_reject_mg); + RUN_TEST(test_projection_zero_init_backward_compat); + + printf("\n========================================\n"); + return UNITY_END(); +} From ea02b828ed6d38bc618b0973af05b60efa8c98e9 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 18 Jul 2026 08:08:38 +0300 Subject: [PATCH 3/9] Propagate non-INVALID multigrid probe failures in projection_init The MG probe-init previously mapped any failure to CFD_ERROR_UNSUPPORTED with a grid-dimension message, masking other causes (e.g. CFD_ERROR_NOMEM) and overwriting their error context. Only the non-2^k+1 grid rejection (CFD_ERROR_INVALID) is remapped now; other statuses propagate unchanged so the real failure cause survives. --- lib/src/api/solver_registry.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/lib/src/api/solver_registry.c b/lib/src/api/solver_registry.c index 2d6d437d..5ad5b1ea 100644 --- a/lib/src/api/solver_registry.c +++ b/lib/src/api/solver_registry.c @@ -939,11 +939,18 @@ static cfd_status_t projection_init(ns_solver_t* solver, const grid* grid, const grid->dx[0], grid->dy[0], (grid->nz > 1 && grid->dz) ? grid->dz[0] : 0.0, NULL); poisson_solver_destroy(probe); - if (probe_status != CFD_SUCCESS) { + if (probe_status == CFD_ERROR_INVALID) { + /* Only the grid-dimension rejection (non-2^k+1) is remapped to + * UNSUPPORTED so the caller can pick a different pressure solver. */ cfd_set_error(CFD_ERROR_UNSUPPORTED, "Multigrid pressure solver requires 2^k+1 grid points per active dimension"); return CFD_ERROR_UNSUPPORTED; } + if (probe_status != CFD_SUCCESS) { + /* Other failures (e.g. CFD_ERROR_NOMEM) propagate unchanged so the + * real cause and its error context are not masked. */ + return probe_status; + } } projection_context* ctx = (projection_context*)cfd_malloc(sizeof(projection_context)); if (!ctx) { From 3d1e0433005b5b16cd671ecf8f0a63db824fe711 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 18 Jul 2026 08:08:56 +0300 Subject: [PATCH 4/9] Make the zero-init pressure-solver test isolate the selection The previous "zero-init backward compat" test ran make_params() twice, so both configs went through ns_solver_params_default() and only differed by NS_PRESSURE_SOLVER_DEFAULT vs (t)0 - the same value - making the memcmp assertions compare identical configurations. Build three configs from one genuinely zero-initialized struct that differ only in pressure_solver: a zero-init run must match the explicit-CG run bitwise and must differ from the explicit-MG run, proving the zero value selects CG (not the multigrid path). Anchoring on a zero base avoids default()'s nonzero source amplitudes, which otherwise confound the compare. The helper now takes a full ns_solver_params_t so callers control exactly how each field was initialized. --- .../cpu/test_projection_pressure_solver.c | 85 ++++++++++++++----- 1 file changed, 62 insertions(+), 23 deletions(-) diff --git a/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c b/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c index aeeb3a8b..29e378ab 100644 --- a/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c +++ b/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c @@ -53,11 +53,13 @@ static ns_solver_params_t make_params(ns_pressure_solver_t pressure_solver) { /** * Run the scalar "projection" solver for NUM_STEPS on a Taylor-Green field. - * Returns the init status; on CFD_SUCCESS the field holds the final state. + * Takes a fully populated params struct so callers can exercise exactly how + * the fields were initialized (default() vs zero-init). Returns the init + * status; on CFD_SUCCESS the field holds the final state. */ -static cfd_status_t run_projection_taylor_green(ns_pressure_solver_t pressure_solver, +static cfd_status_t run_projection_taylor_green(const ns_solver_params_t* params_in, flow_field* field, const grid* g) { - ns_solver_params_t params = make_params(pressure_solver); + ns_solver_params_t params = *params_in; ns_solver_registry_t* registry = cfd_registry_create(); TEST_ASSERT_NOT_NULL_MESSAGE(registry, "Failed to create registry"); @@ -141,12 +143,15 @@ void test_projection_mg_matches_cg(void) { TEST_ASSERT_NOT_NULL(field_mg); TEST_ASSERT_NOT_NULL(field_pcg); + ns_solver_params_t params_cg = make_params(NS_PRESSURE_SOLVER_DEFAULT); + ns_solver_params_t params_mg = make_params(NS_PRESSURE_SOLVER_MULTIGRID); + ns_solver_params_t params_pcg = make_params(NS_PRESSURE_SOLVER_PCG_MG); TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, - run_projection_taylor_green(NS_PRESSURE_SOLVER_DEFAULT, field_cg, g)); + run_projection_taylor_green(¶ms_cg, field_cg, g)); TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, - run_projection_taylor_green(NS_PRESSURE_SOLVER_MULTIGRID, field_mg, g)); + run_projection_taylor_green(¶ms_mg, field_mg, g)); TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, - run_projection_taylor_green(NS_PRESSURE_SOLVER_PCG_MG, field_pcg, g)); + run_projection_taylor_green(¶ms_pcg, field_pcg, g)); size_t n = (size_t)GRID_N * GRID_N; double u_norm = test_compute_l2_norm(field_cg->u, n); @@ -274,7 +279,13 @@ void test_projection_backends_reject_mg(void) { void test_projection_zero_init_backward_compat(void) { printf("\n Testing zero-init params preserve default CG behavior...\n"); + /* Backward-compat contract: a caller that predates the pressure_solver + * field - one that zero-initializes ns_solver_params_t and never touches + * pressure_solver - must still get the original CG pressure solve. That + * only holds because the CG path is selected by enum value 0. */ ns_solver_params_t defaults = ns_solver_params_default(); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, (int)NS_PRESSURE_SOLVER_DEFAULT, + "Zero-init safety requires the default CG path to map to enum value 0"); TEST_ASSERT_EQUAL_INT_MESSAGE(NS_PRESSURE_SOLVER_DEFAULT, defaults.pressure_solver, "ns_solver_params_default must leave pressure_solver at the zero default"); @@ -282,28 +293,56 @@ void test_projection_zero_init_backward_compat(void) { TEST_ASSERT_NOT_NULL(g); grid_initialize_uniform(g); - flow_field* field_a = flow_field_create(GRID_N, GRID_N, 1); - flow_field* field_b = flow_field_create(GRID_N, GRID_N, 1); - TEST_ASSERT_NOT_NULL(field_a); - TEST_ASSERT_NOT_NULL(field_b); + flow_field* field_zero = flow_field_create(GRID_N, GRID_N, 1); + flow_field* field_cg = flow_field_create(GRID_N, GRID_N, 1); + flow_field* field_mg = flow_field_create(GRID_N, GRID_N, 1); + TEST_ASSERT_NOT_NULL(field_zero); + TEST_ASSERT_NOT_NULL(field_cg); + TEST_ASSERT_NOT_NULL(field_mg); - /* Same run twice: once relying on the default, once with the explicit - * enum value. The scalar path is deterministic, so bitwise equality. */ + /* Build three configs from one genuinely zero-initialized struct (the + * legacy caller pattern) that differ ONLY in pressure_solver. Anchoring on + * a common zero base keeps every other field identical, so any difference + * in the result is attributable solely to the pressure-solver selection - + * not to default()'s nonzero source amplitudes or other tuning fields. */ + ns_solver_params_t base; + memset(&base, 0, sizeof(base)); + base.dt = TEST_DT; + base.mu = 0.01; + base.max_iter = 1; + TEST_ASSERT_EQUAL_INT_MESSAGE(NS_PRESSURE_SOLVER_DEFAULT, base.pressure_solver, + "A zero-initialized params struct must leave pressure_solver at the CG default"); + + ns_solver_params_t params_zero = base; /* pressure_solver left at zero-init 0 */ + ns_solver_params_t params_cg = base; + params_cg.pressure_solver = NS_PRESSURE_SOLVER_DEFAULT; /* explicit CG */ + ns_solver_params_t params_mg = base; + params_mg.pressure_solver = NS_PRESSURE_SOLVER_MULTIGRID; /* explicit MG */ + + TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, + run_projection_taylor_green(¶ms_zero, field_zero, g)); TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, - run_projection_taylor_green(NS_PRESSURE_SOLVER_DEFAULT, field_a, g)); + run_projection_taylor_green(¶ms_cg, field_cg, g)); TEST_ASSERT_EQUAL_INT(CFD_SUCCESS, - run_projection_taylor_green((ns_pressure_solver_t)0, field_b, g)); + run_projection_taylor_green(¶ms_mg, field_mg, g)); + /* The zero-init run must take the CG path: bitwise identical to the + * explicit-CG run (scalar solve is deterministic), and NOT the multigrid + * run. The second check guards against a regression that maps the zero + * value to some other pressure solver. */ size_t bytes = (size_t)GRID_N * GRID_N * sizeof(double); - TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_a->u, field_b->u, bytes), - "u fields must be bitwise identical"); - TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_a->v, field_b->v, bytes), - "v fields must be bitwise identical"); - TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_a->p, field_b->p, bytes), - "p fields must be bitwise identical"); - - flow_field_destroy(field_a); - flow_field_destroy(field_b); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_zero->u, field_cg->u, bytes), + "zero-init u field must match explicit CG bitwise"); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_zero->v, field_cg->v, bytes), + "zero-init v field must match explicit CG bitwise"); + TEST_ASSERT_EQUAL_INT_MESSAGE(0, memcmp(field_zero->p, field_cg->p, bytes), + "zero-init p field must match explicit CG bitwise"); + TEST_ASSERT_TRUE_MESSAGE(memcmp(field_zero->u, field_mg->u, bytes) != 0, + "zero-init must select CG, not the multigrid pressure solver"); + + flow_field_destroy(field_zero); + flow_field_destroy(field_cg); + flow_field_destroy(field_mg); grid_destroy(g); } From bc48384c065fa2e6aa6661ee823d9fd97cf6bf76 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 18 Jul 2026 08:09:06 +0300 Subject: [PATCH 5/9] Note that L2=-1.0 is a non-convergence sentinel in the tuning example poisson_solver_tuning.c prints L2=-1.0 when a solve does not return CFD_SUCCESS (max_iter exit), so the Jacobi row's L2=-1.0e+00 is a "not computed" marker rather than a real negative error. Call this out in the Expected Output so readers don't misread it as a valid value. --- docs/guides/examples.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/guides/examples.md b/docs/guides/examples.md index 2fa2584b..203fe766 100644 --- a/docs/guides/examples.md +++ b/docs/guides/examples.md @@ -704,6 +704,11 @@ CFD Platform Diagnostics BiCGSTAB 1 res=6.3e-11 L2=1.0e-04 0 ms converged ``` +> **Note:** `L2=-1.0e+00` is a "not computed" sentinel, not a real error value. The +> example only computes the L2 error when the solve returns `CFD_SUCCESS`; on a +> `max_iter` (non-converged) exit it prints `-1.0` instead. Jacobi does not converge +> on this RHS within 10000 iterations, so its L2 error is reported as `-1.0`. + --- ### 14. poiseuille_stretched_grid.c From faed1e00b1f5748acab39ab06cff4b4483db6883 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 18 Jul 2026 20:29:49 +0300 Subject: [PATCH 6/9] Use the real projection_gpu solver name in the reference docs The reference docs advertised a projection_jacobi_gpu solver type and an NS_SOLVER_TYPE_PROJECTION_JACOBI_GPU macro, neither of which exists: the registry registers NS_SOLVER_TYPE_PROJECTION_GPU ("projection_gpu") and the GPU source is solver_projection_gpu.cu. Copy-pasting the documented name would fail at cfd_solver_create(). Also corrects the backend table description - the GPU projection uses a Conjugate Gradient pressure solve, not Jacobi. Pre-existing: master already carried the stale name; this PR made it visible by adding adjacent text that uses the correct one. Historical records under technical-notes/ and validation/ are left as-is. --- docs/reference/api-reference.md | 2 +- docs/reference/solvers.md | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/reference/api-reference.md b/docs/reference/api-reference.md index 17a2a4e4..44298846 100644 --- a/docs/reference/api-reference.md +++ b/docs/reference/api-reference.md @@ -756,7 +756,7 @@ cfd_aligned_free(data); // Correct #define NS_SOLVER_TYPE_PROJECTION_OPTIMIZED "projection_optimized" #define NS_SOLVER_TYPE_PROJECTION_OMP "projection_omp" #define NS_SOLVER_TYPE_RK2 "rk2" -#define NS_SOLVER_TYPE_PROJECTION_JACOBI_GPU "projection_jacobi_gpu" +#define NS_SOLVER_TYPE_PROJECTION_GPU "projection_gpu" ``` ## Version Information diff --git a/docs/reference/solvers.md b/docs/reference/solvers.md index 6e4978a9..e60e8e1e 100644 --- a/docs/reference/solvers.md +++ b/docs/reference/solvers.md @@ -76,7 +76,7 @@ Chorin's projection method - properly enforces incompressibility constraint. | `projection` | Scalar | Basic implementation | | `projection_optimized` | SIMD | SIMD-optimized (runtime detection: AVX2/NEON) | | `projection_omp` | OpenMP | Multi-threaded | -| `projection_jacobi_gpu` | GPU | CUDA-accelerated (Jacobi iteration) | +| `projection_gpu` | GPU | CUDA-accelerated (CG pressure solve) | **Pressure solver selection** (`ns_solver_params_t.pressure_solver`): @@ -463,7 +463,7 @@ GPU-accelerated using CUDA: // Check if GPU should be used gpu_config_t config = gpu_config_default(); if (gpu_should_use(&config, nx, ny, num_steps)) { - solver = cfd_solver_create(registry, "projection_jacobi_gpu"); + solver = cfd_solver_create(registry, "projection_gpu"); } else { solver = cfd_solver_create(registry, "projection_optimized"); } @@ -481,7 +481,7 @@ if (gpu_should_use(&config, nx, ny, num_steps)) { | projection | 19.0 | 1.0x | High | | projection_optimized | 5.3 | 3.6x | High | | projection_omp (8 cores) | 4.2 | 4.5x | High | -| projection_jacobi_gpu | 8.4 | 0.45x† | High | +| projection_gpu | 8.4 | 0.45x† | High | † GPU slower on small grids due to data transfer overhead @@ -499,7 +499,7 @@ if (gpu_should_use(&config, nx, ny, num_steps)) { | Solver | Time (s) | Speedup | |--------|----------|---------| | projection_optimized | 824 | 1.0x | -| projection_jacobi_gpu | 68 | 12.1x | +| projection_gpu | 68 | 12.1x | ## Choosing a Solver @@ -515,7 +515,7 @@ Need strict incompressibility enforcement? └─ Yes → Use Projection Method family ├─ Small grid (<100×100) → projection ├─ Medium grid (100-500) → projection_optimized or projection_omp - └─ Large grid (>500) → projection_jacobi_gpu + └─ Large grid (>500) → projection_gpu GPU available and grid >200×200? └─ Use CUDA variant for 10-50x speedup @@ -530,7 +530,7 @@ GPU available and grid >200×200? **Production Simulations:** - `projection_optimized` or `projection_omp` (medium grids) -- `projection_jacobi_gpu` (large grids) +- `projection_gpu` (large grids) - Best accuracy and performance **Benchmarking:** From bc6f985a64a4a1fe6abbb94d716463a422f56437 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 18 Jul 2026 20:34:44 +0300 Subject: [PATCH 7/9] Set an error message when the multigrid probe factory is unavailable poisson_solver_create() returning NULL made projection_init return CFD_ERROR_UNSUPPORTED with no last-error message, leaving callers and tests without a reason for the rejection. Set an explicit message, matching the other rejection paths in this init. --- lib/src/api/solver_registry.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/src/api/solver_registry.c b/lib/src/api/solver_registry.c index 5ad5b1ea..b0671720 100644 --- a/lib/src/api/solver_registry.c +++ b/lib/src/api/solver_registry.c @@ -932,6 +932,8 @@ static cfd_status_t projection_init(ns_solver_t* solver, const grid* grid, const poisson_solver_t* probe = poisson_solver_create( POISSON_METHOD_MULTIGRID, POISSON_BACKEND_SCALAR); if (!probe) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, + "Multigrid Poisson solver is unavailable in this build"); return CFD_ERROR_UNSUPPORTED; } cfd_status_t probe_status = poisson_solver_init( From b6703a5f668dee8e17c48e95901d362bbac7f174 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 18 Jul 2026 22:12:38 +0300 Subject: [PATCH 8/9] Stop the multigrid probe from allocating the full hierarchy projection_init probe-initializes a multigrid solver purely to validate the grid dimensions, but a default init builds the entire coarse-grid hierarchy on conforming grids - an O(N) allocation and transient memory spike that the real pressure solve then rebuilds from scratch. Pass mg_max_levels = 1 so the probe stops at the finest level. Level 0 borrows the caller's buffers and needs no residual buffer, so the probe now allocates only the context and one level descriptor. The 2^k+1 dimension check runs before any of that, so the rejection behavior is unchanged. --- lib/src/api/solver_registry.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/lib/src/api/solver_registry.c b/lib/src/api/solver_registry.c index b0671720..7edda312 100644 --- a/lib/src/api/solver_registry.c +++ b/lib/src/api/solver_registry.c @@ -936,10 +936,17 @@ static cfd_status_t projection_init(ns_solver_t* solver, const grid* grid, const "Multigrid Poisson solver is unavailable in this build"); return CFD_ERROR_UNSUPPORTED; } + /* The probe only needs the dimension checks, which multigrid init runs + * before it builds anything. mg_max_levels = 1 stops it from + * allocating the coarse-grid hierarchy that the real pressure solve + * would rebuild anyway, avoiding an O(N) transient spike at init. */ + poisson_solver_params_t probe_params = poisson_solver_params_default(); + probe_params.mg_max_levels = 1; + cfd_status_t probe_status = poisson_solver_init( probe, grid->nx, grid->ny, grid->nz, grid->dx[0], grid->dy[0], - (grid->nz > 1 && grid->dz) ? grid->dz[0] : 0.0, NULL); + (grid->nz > 1 && grid->dz) ? grid->dz[0] : 0.0, &probe_params); poisson_solver_destroy(probe); if (probe_status == CFD_ERROR_INVALID) { /* Only the grid-dimension rejection (non-2^k+1) is remapped to From 5a3fa787feec05c3d4168d2d6f347adcb82a8958 Mon Sep 17 00:00:00 2001 From: shaia Date: Fri, 24 Jul 2026 20:44:25 +0300 Subject: [PATCH 9/9] Keep degenerate grids as INVALID, not UNSUPPORTED, for MG pressure solve The INVALID->UNSUPPORTED remap in projection_init assumed the probe's CFD_ERROR_INVALID could only mean a non-2^k+1 grid. But poisson_solver_init also returns INVALID for degenerate grids (nx<3, ny<3, or nz==2), so those were misreported as UNSUPPORTED with a "requires 2^k+1" message. Screen degenerate grids up front and return INVALID with a fitting message. The two statuses drive different recovery: UNSUPPORTED tells the caller to try another pressure solver, which is useless advice for a grid too small for any solver. After this screen the probe's INVALID is unambiguous. Adds test_projection_mg_rejects_degenerate_grid_as_invalid covering a 2x2 grid, alongside the existing non-2^k+1 (UNSUPPORTED) case. --- lib/src/api/solver_registry.c | 12 +++++ .../cpu/test_projection_pressure_solver.c | 53 ++++++++++++++++++- 2 files changed, 63 insertions(+), 2 deletions(-) diff --git a/lib/src/api/solver_registry.c b/lib/src/api/solver_registry.c index 7edda312..3a64fd43 100644 --- a/lib/src/api/solver_registry.c +++ b/lib/src/api/solver_registry.c @@ -927,6 +927,18 @@ static cfd_status_t projection_init(ns_solver_t* solver, const grid* grid, const params->pressure_solver != NS_PRESSURE_SOLVER_PCG_MG) { return CFD_ERROR_INVALID; } + /* Screen degenerate grids first. poisson_solver_init reports those + * with CFD_ERROR_INVALID as well, so rejecting them here keeps the + * probe's INVALID unambiguous: below, it can only mean the + * multigrid-specific non-2^k+1 rejection. A grid too small to hold an + * interior cell is a caller error, not an unsupported configuration, + * so it stays INVALID rather than being remapped to UNSUPPORTED. */ + if (grid->nx < 3 || grid->ny < 3 || (grid->nz > 1 && grid->nz < 3)) { + cfd_set_error(CFD_ERROR_INVALID, + "Multigrid pressure solver requires at least 3 points per active dimension"); + return CFD_ERROR_INVALID; + } + /* Both MG modes require a multigrid hierarchy on this exact grid: * probe-init and reject non-2^k+1 dimensions up front. */ poisson_solver_t* probe = poisson_solver_create( diff --git a/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c b/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c index 29e378ab..b8ddec55 100644 --- a/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c +++ b/tests/solvers/navier_stokes/cpu/test_projection_pressure_solver.c @@ -7,6 +7,8 @@ * - The MG modes agree with the default CG pressure solve within solver * tolerance * - Non-2^k+1 grids are rejected at init with CFD_ERROR_UNSUPPORTED + * - Degenerate grids are rejected with CFD_ERROR_INVALID instead, since no + * pressure solver choice can rescue them * - projection_optimized / projection_omp / projection_gpu reject any * non-default selection with CFD_ERROR_UNSUPPORTED (no silent fallbacks) * - Zero-initialized params keep the existing CG behavior bit-for-bit @@ -37,8 +39,9 @@ void tearDown(void) { cfd_finalize(); } -#define GRID_N 33 /* 2^5+1: multigrid-conforming */ -#define GRID_N_BAD 30 /* not 2^k+1 */ +#define GRID_N 33 /* 2^5+1: multigrid-conforming */ +#define GRID_N_BAD 30 /* not 2^k+1, but large enough to be a valid grid */ +#define GRID_N_DEGENERATE 2 /* too small to hold an interior cell */ #define NUM_STEPS 5 #define TEST_DT 1e-3 @@ -185,6 +188,51 @@ void test_projection_mg_matches_cg(void) { grid_destroy(g); } +//============================================================================= +// TEST: DEGENERATE GRID IS REJECTED AS INVALID, NOT UNSUPPORTED +//============================================================================= + +/** + * A grid too small to hold an interior cell is a caller error (INVALID), not + * a configuration this solver merely lacks support for (UNSUPPORTED). The two + * statuses drive different recovery: UNSUPPORTED tells the caller to pick a + * different pressure solver, which would be useless advice here since no + * pressure solver can work on a 2x2 grid. + */ +void test_projection_mg_rejects_degenerate_grid_as_invalid(void) { + printf("\n Testing MG pressure solver rejects degenerate grids as INVALID...\n"); + + grid* g = grid_create(GRID_N_DEGENERATE, GRID_N_DEGENERATE, 1, + 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + ns_solver_registry_t* registry = cfd_registry_create(); + TEST_ASSERT_NOT_NULL(registry); + cfd_registry_register_defaults(registry); + + ns_pressure_solver_t modes[] = { + NS_PRESSURE_SOLVER_MULTIGRID, + NS_PRESSURE_SOLVER_PCG_MG, + }; + for (size_t m = 0; m < sizeof(modes) / sizeof(modes[0]); m++) { + ns_solver_t* slv = cfd_solver_create(registry, NS_SOLVER_TYPE_PROJECTION); + TEST_ASSERT_NOT_NULL_MESSAGE(slv, "projection solver not available"); + + ns_solver_params_t params = make_params(modes[m]); + cfd_status_t init_status = solver_init(slv, g, ¶ms); + printf(" mode %d on %dx%d: init status %d\n", + (int)modes[m], GRID_N_DEGENERATE, GRID_N_DEGENERATE, (int)init_status); + TEST_ASSERT_EQUAL_INT_MESSAGE(CFD_ERROR_INVALID, init_status, + "MG pressure solver on a degenerate grid must fail init with INVALID"); + + solver_destroy(slv); + } + + cfd_registry_destroy(registry); + grid_destroy(g); +} + //============================================================================= // TEST: NON-2^K+1 GRID IS REJECTED AT INIT //============================================================================= @@ -360,6 +408,7 @@ int main(void) { RUN_TEST(test_projection_mg_divergence_free); RUN_TEST(test_projection_pcg_mg_divergence_free); RUN_TEST(test_projection_mg_matches_cg); + RUN_TEST(test_projection_mg_rejects_degenerate_grid_as_invalid); RUN_TEST(test_projection_mg_rejects_non_pow2_grid); RUN_TEST(test_projection_backends_reject_mg); RUN_TEST(test_projection_zero_init_backward_compat);