From 03c43456b36c7c98c1ca3a50ec7883a11d8413b3 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 25 Jul 2026 16:21:38 +0300 Subject: [PATCH 1/6] Add Jacobi and BiCGSTAB OpenMP linear solvers Completes the OpenMP linear-solver tier (previously only CG, GMRES, and Red-Black SOR), so the parallel projection/pressure paths can select these methods instead of falling back to unavailable. Both are selectable via poisson_solver_create(method, POISSON_BACKEND_OMP) and follow the existing cg_omp/redblack_omp pattern: per-element updates are identical to the scalar reference, so results agree within reduction rounding (verified by new scalar-vs-OMP consistency tests). BiCGSTAB keeps the scalar solve loop and breakdown checks verbatim, parallelizing only the primitives. Plain lexicographic SOR stays scalar/SIMD-only: it is inherently sequential, and its parallel form is the existing Red-Black SOR OMP solver. Also updates the stale BiCGSTAB backend test, which asserted OMP create() returns NULL, to the runtime availability contract. --- CHANGELOG.md | 9 + ROADMAP.md | 14 +- lib/CMakeLists.txt | 2 + lib/include/cfd/solvers/poisson_solver.h | 2 + lib/src/solvers/linear/linear_solver.c | 8 + .../solvers/linear/linear_solver_internal.h | 2 + .../linear/omp/linear_solver_bicgstab_omp.c | 592 ++++++++++++++++++ .../linear/omp/linear_solver_jacobi_omp.c | 160 +++++ tests/math/test_bicgstab.c | 23 +- tests/math/test_omp_consistency.c | 199 ++++++ 10 files changed, 1000 insertions(+), 11 deletions(-) create mode 100644 lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c create mode 100644 lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c diff --git a/CHANGELOG.md b/CHANGELOG.md index ad0cd3bc..0dae69d3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- **Jacobi and BiCGSTAB OpenMP backends** — completes the OpenMP linear-solver + tier (previously CG, GMRES, and Red-Black SOR). Both are selectable via + `poisson_solver_create(method, POISSON_BACKEND_OMP)`; per-element updates are + identical to the scalar reference, so results match within reduction rounding + (verified by scalar-vs-OMP consistency tests). Plain lexicographic SOR remains + scalar/SIMD-only — its parallel form is the existing Red-Black SOR OMP solver + (`lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c`, + `lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c`, + `tests/math/test_omp_consistency.c`). - **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 diff --git a/ROADMAP.md b/ROADMAP.md index 88a57a38..50fcc003 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -48,15 +48,20 @@ The single source of truth for backend gaps. Each algorithm targets scalar (CPU) | | RK4 (classical)| done | done | — | done | done | | **Energy Eq.** | Advec-diff + Boussinesq + thermal BCs | done | done | — | done | done | | **Turbulence** | k-ε / SA + wall functions | done | done | — | done | — | -| **Linear Solvers** | Jacobi | done | done | done | — | done | -| | SOR | done | done | done | — | done | +| **Linear Solvers** | Jacobi | done | done | done | done | done | +| | SOR | done | done | done | —¹ | done | | | Red-Black SOR | done | done | done | done | done | | | CG / PCG | done | done | done | done | done | -| | BiCGSTAB | done | done | done | — | done | +| | BiCGSTAB | done | done | done | done | done | | | GMRES(m) | done | done | done | done | — | | | Multigrid (GMG)| done | — | — | — | — | | **Boundary Conds** | All types | done | done | done | done | done | +¹ Plain (lexicographic) SOR is inherently sequential — each update reads +already-updated neighbors. Its parallel form is **Red-Black SOR**, which has an +OMP backend; a "plain SOR OMP" would either change the numerics silently or need +low-value wavefront machinery, so it is intentionally omitted. + ### Known Limitations Genuine constraints to be aware of (not backlog items): @@ -127,6 +132,9 @@ non-symmetric operators arrive, e.g. implicit advection-diffusion in §1.5). - [x] GMRES (Generalized Minimal Residual) for non-symmetric systems — scalar, AVX2, NEON, OMP (done). GPU variant deferred. - [ ] GMRES GPU backend (deferred from the initial GMRES landing) +- [x] Jacobi and BiCGSTAB OpenMP backends — completes the OMP linear tier + (Jacobi, Red-Black SOR, CG/PCG, BiCGSTAB, GMRES). Plain lexicographic SOR + stays scalar/SIMD-only; its parallel form is Red-Black SOR OMP. - [ ] SSOR (Symmetric SOR) preconditioner - [ ] ILU preconditioner - [x] Geometric multigrid — scalar backend; V/W/F(FMG) cycles, Red-Black GS or weighted diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index b0ea9557..d1d66d61 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -185,8 +185,10 @@ set(CFD_OMP_SOURCES # Turbulence models (RANS) - OMP src/solvers/turbulence/omp/turbulence_solver_omp.c # Linear solvers - OMP + src/solvers/linear/omp/linear_solver_jacobi_omp.c src/solvers/linear/omp/linear_solver_redblack_omp.c src/solvers/linear/omp/linear_solver_cg_omp.c + src/solvers/linear/omp/linear_solver_bicgstab_omp.c src/solvers/linear/omp/linear_solver_gmres_omp.c ) diff --git a/lib/include/cfd/solvers/poisson_solver.h b/lib/include/cfd/solvers/poisson_solver.h index 1887ab62..6396bf66 100644 --- a/lib/include/cfd/solvers/poisson_solver.h +++ b/lib/include/cfd/solvers/poisson_solver.h @@ -426,6 +426,7 @@ CFD_LIBRARY_EXPORT bool poisson_solver_backend_available(poisson_solver_backend_ * ============================================================================ */ #define POISSON_SOLVER_TYPE_JACOBI_SCALAR "jacobi_scalar" +#define POISSON_SOLVER_TYPE_JACOBI_OMP "jacobi_omp" #define POISSON_SOLVER_TYPE_JACOBI_SIMD "jacobi_simd" #define POISSON_SOLVER_TYPE_SOR_SCALAR "sor_scalar" #define POISSON_SOLVER_TYPE_SOR_SIMD "sor_simd" @@ -436,6 +437,7 @@ CFD_LIBRARY_EXPORT bool poisson_solver_backend_available(poisson_solver_backend_ #define POISSON_SOLVER_TYPE_CG_OMP "cg_omp" #define POISSON_SOLVER_TYPE_CG_SIMD "cg_simd" #define POISSON_SOLVER_TYPE_BICGSTAB_SCALAR "bicgstab_scalar" +#define POISSON_SOLVER_TYPE_BICGSTAB_OMP "bicgstab_omp" #define POISSON_SOLVER_TYPE_BICGSTAB_SIMD "bicgstab_simd" #define POISSON_SOLVER_TYPE_BICGSTAB_GPU "bicgstab_gpu" #define POISSON_SOLVER_TYPE_GMRES_SCALAR "gmres_scalar" diff --git a/lib/src/solvers/linear/linear_solver.c b/lib/src/solvers/linear/linear_solver.c index 1c845743..bae17ade 100644 --- a/lib/src/solvers/linear/linear_solver.c +++ b/lib/src/solvers/linear/linear_solver.c @@ -181,6 +181,10 @@ poisson_solver_t* poisson_solver_create( switch (backend) { case POISSON_BACKEND_SIMD: return create_jacobi_simd_solver(); +#ifdef CFD_ENABLE_OPENMP + case POISSON_BACKEND_OMP: + return create_jacobi_omp_solver(); +#endif #ifdef CFD_HAS_CUDA case POISSON_BACKEND_GPU: return create_jacobi_gpu_solver(); @@ -246,6 +250,10 @@ poisson_solver_t* poisson_solver_create( switch (backend) { case POISSON_BACKEND_SIMD: return create_bicgstab_simd_solver(); +#ifdef CFD_ENABLE_OPENMP + case POISSON_BACKEND_OMP: + return create_bicgstab_omp_solver(); +#endif #ifdef CFD_HAS_CUDA case POISSON_BACKEND_GPU: return create_bicgstab_gpu_solver(); diff --git a/lib/src/solvers/linear/linear_solver_internal.h b/lib/src/solvers/linear/linear_solver_internal.h index 9d003f7d..b9fc3a6a 100644 --- a/lib/src/solvers/linear/linear_solver_internal.h +++ b/lib/src/solvers/linear/linear_solver_internal.h @@ -41,8 +41,10 @@ poisson_solver_t* create_redblack_scalar_solver(void); poisson_solver_t* create_redblack_simd_solver(void); #ifdef CFD_ENABLE_OPENMP +poisson_solver_t* create_jacobi_omp_solver(void); poisson_solver_t* create_redblack_omp_solver(void); poisson_solver_t* create_cg_omp_solver(void); +poisson_solver_t* create_bicgstab_omp_solver(void); #endif /* Conjugate Gradient solvers */ diff --git a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c new file mode 100644 index 00000000..b4657638 --- /dev/null +++ b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c @@ -0,0 +1,592 @@ +/** + * @file linear_solver_bicgstab_omp.c + * @brief BiCGSTAB solver - OpenMP parallelized implementation + * + * Same algorithm as scalar BiCGSTAB (van der Vorst, 1992) but with the vector + * primitives (dot products, axpy, Laplacian, and the per-iteration vector + * updates) parallelized across OpenMP threads. The high-level solve loop and + * all breakdown checks are identical to the scalar reference — only the + * primitives differ, per the algorithm/primitive separation rule. + * + * Each interior element is computed independently and identically regardless of + * thread count; only the dot-product reductions accumulate in a different order, + * so OMP results agree with scalar to within reduction rounding. + */ + +#include "../linear_solver_internal.h" + +#include "cfd/core/indexing.h" +#include "cfd/core/logging.h" +#include "cfd/core/memory.h" + +#include + +#ifdef CFD_ENABLE_OPENMP + +#include + +/* ============================================================================ + * BICGSTAB CONTEXT + * ============================================================================ */ + +typedef struct { + double dx2; /* dx^2 */ + double dy2; /* dy^2 */ + double inv_dz2; /* 1/dz^2 (0 for 2D) */ + + size_t stride_z; /* nx*ny for 3D, 0 for 2D */ + size_t k_start; /* first interior k index */ + size_t k_end; /* one-past-last interior k index */ + + /* BiCGSTAB working vectors (allocated during init) */ + double* r; /* Residual vector */ + double* r_hat; /* Shadow residual (typically r_0) */ + double* p; /* Search direction */ + double* v; /* A * p */ + double* s; /* Intermediate residual */ + double* t; /* A * s */ + + int initialized; +} bicgstab_omp_context_t; + +/* ============================================================================ + * OMP-PARALLELIZED PRIMITIVES + * ============================================================================ */ + +static double dot_product_omp(const double* a, const double* b, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + double sum = 0.0; + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) reduction(+:sum) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + sum += a[idx] * b[idx]; + } + } + } + return sum; +} + +/* y = y + alpha * x (interior points only) */ +static void axpy_omp(double alpha, const double* x, double* y, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + y[idx] += alpha * x[idx]; + } + } + } +} + +/* + * Apply negative Laplacian operator: Ap = -nabla^2(p) + * For Poisson equation nabla^2(x) = rhs we solve -nabla^2(x) = -rhs, which is + * SPD with positive eigenvalues. + */ +static void apply_laplacian_omp(const double* p, double* Ap, + size_t nx, size_t ny, + double dx2, double dy2, double inv_dz2, + size_t k_start, size_t k_end, size_t stride_z) { + double dx2_inv = 1.0 / dx2; + double dy2_inv = 1.0 / dy2; + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + double laplacian = + (p[idx + 1] - 2.0 * p[idx] + p[idx - 1]) * dx2_inv + + (p[idx + nx] - 2.0 * p[idx] + p[idx - nx]) * dy2_inv + + (p[idx + stride_z] + p[idx - stride_z] - 2.0 * p[idx]) * inv_dz2; + Ap[idx] = -laplacian; + } + } + } +} + +/* + * Compute initial residual: r = b - A*x + * For our formulation: r = -rhs - (-nabla^2 x) = -rhs + nabla^2 x + */ +static void compute_residual_omp(const double* x, const double* rhs, double* r, + size_t nx, size_t ny, + double dx2, double dy2, double inv_dz2, + size_t k_start, size_t k_end, size_t stride_z) { + double dx2_inv = 1.0 / dx2; + double dy2_inv = 1.0 / dy2; + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + double laplacian = + (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) * dx2_inv + + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) * dy2_inv + + (x[idx + stride_z] + x[idx - stride_z] - 2.0 * x[idx]) * inv_dz2; + r[idx] = -rhs[idx] + laplacian; + } + } + } +} + +/* dst = src (interior points only) */ +static void copy_vector_omp(const double* src, double* dst, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + dst[idx] = src[idx]; + } + } + } +} + +/* v = 0 (interior points only) */ +static void zero_vector_omp(double* v, size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + v[idx] = 0.0; + } + } + } +} + +/* + * Fused per-iteration vector updates. Each preserves the scalar reference's + * exact per-element floating-point operation order (only the dot-product + * reductions above differ across threads), so OMP and scalar agree tightly. + */ + +/* p = r + beta * (p - omega * v) */ +static void update_p_omp(double* p, const double* r, const double* v, + double beta, double omega, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + p[idx] = r[idx] + beta * (p[idx] - omega * v[idx]); + } + } + } +} + +/* s = r - alpha * v */ +static void update_s_omp(double* s, const double* r, const double* v, + double alpha, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + s[idx] = r[idx] - alpha * v[idx]; + } + } + } +} + +/* x = x + alpha * p + omega * s */ +static void update_x_omp(double* x, const double* p, const double* s, + double alpha, double omega, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + x[idx] += alpha * p[idx] + omega * s[idx]; + } + } + } +} + +/* r = s - omega * t */ +static void update_r_omp(double* r, const double* s, const double* t, + double omega, + size_t nx, size_t ny, + size_t k_start, size_t k_end, size_t stride_z) { + int ny_int = bicgstab_size_to_int(ny); + int nx_int = bicgstab_size_to_int(nx); + + for (size_t k = k_start; k < k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < ny_int - 1; j++) { + for (int i = 1; i < nx_int - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + r[idx] = s[idx] - omega * t[idx]; + } + } + } +} + +/* ============================================================================ + * BICGSTAB OMP IMPLEMENTATION + * ============================================================================ */ + +static cfd_status_t bicgstab_omp_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; + + bicgstab_omp_context_t* ctx = (bicgstab_omp_context_t*)cfd_calloc(1, sizeof(bicgstab_omp_context_t)); + if (!ctx) { + return CFD_ERROR_NOMEM; + } + + ctx->dx2 = dx * dx; + ctx->dy2 = dy * dy; + ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); + poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end); + + /* Allocate working vectors */ + size_t n = nx * ny * nz; + ctx->r = (double*)cfd_calloc(n, sizeof(double)); + ctx->r_hat = (double*)cfd_calloc(n, sizeof(double)); + ctx->p = (double*)cfd_calloc(n, sizeof(double)); + ctx->v = (double*)cfd_calloc(n, sizeof(double)); + ctx->s = (double*)cfd_calloc(n, sizeof(double)); + ctx->t = (double*)cfd_calloc(n, sizeof(double)); + + if (!ctx->r || !ctx->r_hat || !ctx->p || !ctx->v || !ctx->s || !ctx->t) { + cfd_free(ctx->r); + cfd_free(ctx->r_hat); + cfd_free(ctx->p); + cfd_free(ctx->v); + cfd_free(ctx->s); + cfd_free(ctx->t); + cfd_free(ctx); + return CFD_ERROR_NOMEM; + } + + ctx->initialized = 1; + solver->context = ctx; + return CFD_SUCCESS; +} + +static void bicgstab_omp_destroy(poisson_solver_t* solver) { + if (solver && solver->context) { + bicgstab_omp_context_t* ctx = (bicgstab_omp_context_t*)solver->context; + cfd_free(ctx->r); + cfd_free(ctx->r_hat); + cfd_free(ctx->p); + cfd_free(ctx->v); + cfd_free(ctx->s); + cfd_free(ctx->t); + cfd_free(ctx); + solver->context = NULL; + } +} + +/* + * BiCGSTAB solve loop — structurally identical to the scalar reference. It uses + * its own solve function because it maintains complex state across iterations + * (rho, omega, alpha), has two convergence checks per iteration, and requires + * two matrix-vector products per iteration. + */ +static cfd_status_t bicgstab_omp_solve( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + poisson_solver_stats_t* stats) +{ + (void)x_temp; /* BiCGSTAB doesn't use the temp buffer */ + + bicgstab_omp_context_t* ctx = (bicgstab_omp_context_t*)solver->context; + size_t nx = solver->nx; + size_t ny = solver->ny; + double dx2 = ctx->dx2; + double dy2 = ctx->dy2; + double inv_dz2 = ctx->inv_dz2; + size_t stride_z = ctx->stride_z; + size_t k_start = ctx->k_start; + size_t k_end = ctx->k_end; + + double* r = ctx->r; + double* r_hat = ctx->r_hat; + double* p = ctx->p; + double* v = ctx->v; + double* s = ctx->s; + double* t = ctx->t; + + poisson_solver_params_t* params = &solver->params; + double start_time = poisson_solver_get_time_ms(); + + /* Apply initial boundary conditions */ + poisson_solver_apply_bc(solver, x); + + /* Compute initial residual: r_0 = b - A*x_0 */ + compute_residual_omp(x, rhs, r, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + + /* r_hat = r_0 (shadow residual) */ + copy_vector_omp(r, r_hat, nx, ny, k_start, k_end, stride_z); + + /* Initialize: rho = alpha = omega = 1, v = p = 0 */ + double rho = 1.0; + double alpha = 1.0; + double omega = 1.0; + zero_vector_omp(v, nx, ny, k_start, k_end, stride_z); + zero_vector_omp(p, nx, ny, k_start, k_end, stride_z); + + /* Compute initial residual norm */ + double r_dot_r = dot_product_omp(r, r, nx, ny, k_start, k_end, stride_z); + double initial_res = sqrt(r_dot_r); + + if (stats) { + stats->initial_residual = initial_res; + } + + /* Check if already converged */ + double tolerance = params->tolerance * initial_res; + if (tolerance < params->absolute_tolerance) { + tolerance = params->absolute_tolerance; + } + + if (initial_res < params->absolute_tolerance) { + if (stats) { + stats->status = POISSON_CONVERGED; + stats->iterations = 0; + stats->final_residual = initial_res; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_SUCCESS; + } + + int converged = 0; + int iter; + double res_norm = initial_res; + + for (iter = 0; iter < params->max_iterations; iter++) { + /* rho_new = (r_hat, r) */ + double rho_new = dot_product_omp(r_hat, r, nx, ny, k_start, k_end, stride_z); + + /* Check for breakdown */ + if (fabs(rho_new) < BICGSTAB_BREAKDOWN_THRESHOLD) { + if (stats) { + stats->status = POISSON_STAGNATED; + stats->iterations = iter + 1; + stats->final_residual = res_norm; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_ERROR_MAX_ITER; + } + + /* beta = (rho_new / rho) * (alpha / omega) */ + double beta = (rho_new / rho) * (alpha / omega); + + /* p = r + beta * (p - omega * v) */ + update_p_omp(p, r, v, beta, omega, nx, ny, k_start, k_end, stride_z); + + /* v = A * p */ + apply_laplacian_omp(p, v, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + + /* alpha = rho_new / (r_hat, v) */ + double r_hat_dot_v = dot_product_omp(r_hat, v, nx, ny, k_start, k_end, stride_z); + + /* Check for breakdown */ + if (fabs(r_hat_dot_v) < BICGSTAB_BREAKDOWN_THRESHOLD) { + if (stats) { + stats->status = POISSON_STAGNATED; + stats->iterations = iter + 1; + stats->final_residual = res_norm; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_ERROR_MAX_ITER; + } + + alpha = rho_new / r_hat_dot_v; + + /* s = r - alpha * v */ + update_s_omp(s, r, v, alpha, nx, ny, k_start, k_end, stride_z); + + /* Check for early convergence on s */ + double s_norm = sqrt(dot_product_omp(s, s, nx, ny, k_start, k_end, stride_z)); + if (s_norm < tolerance || s_norm < params->absolute_tolerance) { + /* Update x and return */ + axpy_omp(alpha, p, x, nx, ny, k_start, k_end, stride_z); + res_norm = s_norm; + converged = 1; + break; + } + + /* t = A * s */ + apply_laplacian_omp(s, t, nx, ny, dx2, dy2, inv_dz2, k_start, k_end, stride_z); + + /* omega = (t, s) / (t, t) */ + double t_dot_s = dot_product_omp(t, s, nx, ny, k_start, k_end, stride_z); + double t_dot_t = dot_product_omp(t, t, nx, ny, k_start, k_end, stride_z); + + /* Check for breakdown */ + if (fabs(t_dot_t) < BICGSTAB_BREAKDOWN_THRESHOLD) { + /* Update x with available progress */ + axpy_omp(alpha, p, x, nx, ny, k_start, k_end, stride_z); + if (stats) { + stats->status = POISSON_STAGNATED; + stats->iterations = iter + 1; + stats->final_residual = s_norm; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_ERROR_MAX_ITER; + } + + omega = t_dot_s / t_dot_t; + + /* x = x + alpha * p + omega * s */ + update_x_omp(x, p, s, alpha, omega, nx, ny, k_start, k_end, stride_z); + + /* r = s - omega * t */ + update_r_omp(r, s, t, omega, nx, ny, k_start, k_end, stride_z); + + /* Update rho for next iteration */ + rho = rho_new; + + /* Compute residual norm */ + res_norm = sqrt(dot_product_omp(r, r, nx, ny, k_start, k_end, stride_z)); + + /* Check convergence at intervals */ + if (iter % params->check_interval == 0) { + if (params->verbose) { + CFD_LOG_DEBUG("poisson", "BiCGSTAB-OMP Iter %d: residual = %.6e", iter, res_norm); + } + + if (res_norm < tolerance || res_norm < params->absolute_tolerance) { + converged = 1; + break; + } + } + + /* Check for omega breakdown (would cause division by zero next iteration) */ + if (fabs(omega) < BICGSTAB_BREAKDOWN_THRESHOLD) { + if (stats) { + stats->status = POISSON_STAGNATED; + stats->iterations = iter + 1; + stats->final_residual = res_norm; + stats->elapsed_time_ms = poisson_solver_get_time_ms() - start_time; + } + return CFD_ERROR_MAX_ITER; + } + } + + /* Final convergence check */ + if (!converged && (res_norm < tolerance || res_norm < params->absolute_tolerance)) { + converged = 1; + } + + /* Apply final boundary conditions */ + poisson_solver_apply_bc(solver, x); + + double end_time = poisson_solver_get_time_ms(); + + if (stats) { + stats->iterations = (iter < params->max_iterations) ? (iter + 1) : iter; + stats->final_residual = res_norm; + stats->elapsed_time_ms = end_time - start_time; + stats->status = converged ? POISSON_CONVERGED : POISSON_MAX_ITER; + } + + return converged ? CFD_SUCCESS : CFD_ERROR_MAX_ITER; +} + +/* + * Single iteration is not well-defined for BiCGSTAB (it maintains internal + * state). Provide a minimal implementation that returns the current residual. + */ +static cfd_status_t bicgstab_omp_iterate( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + double* residual) +{ + (void)x_temp; + + if (residual) { + *residual = poisson_solver_compute_residual(solver, x, rhs); + } + return CFD_SUCCESS; +} + +/* ============================================================================ + * FACTORY FUNCTION + * ============================================================================ */ + +poisson_solver_t* create_bicgstab_omp_solver(void) { + poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t)); + if (!solver) { + return NULL; + } + + solver->name = POISSON_SOLVER_TYPE_BICGSTAB_OMP; + solver->description = "BiCGSTAB (OpenMP)"; + solver->method = POISSON_METHOD_BICGSTAB; + solver->backend = POISSON_BACKEND_OMP; + solver->params = poisson_solver_params_default(); + + solver->init = bicgstab_omp_init; + solver->destroy = bicgstab_omp_destroy; + solver->solve = bicgstab_omp_solve; /* BiCGSTAB uses custom solve loop */ + solver->iterate = bicgstab_omp_iterate; + solver->apply_bc = NULL; /* Use default Neumann */ + + return solver; +} + +#endif /* CFD_ENABLE_OPENMP */ diff --git a/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c b/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c new file mode 100644 index 00000000..6826ed84 --- /dev/null +++ b/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c @@ -0,0 +1,160 @@ +/** + * @file linear_solver_jacobi_omp.c + * @brief Jacobi iteration solver - OpenMP parallelized implementation + * + * Same algorithm as scalar Jacobi (double-buffered: reads p_old, writes p_new), + * with the interior update loop parallelized across OpenMP threads. Jacobi is + * embarrassingly parallel — each interior update reads only the previous + * iterate, so there is no cross-cell dependency within a sweep. + */ + +#include "../linear_solver_internal.h" + +#include "cfd/core/indexing.h" +#include "cfd/core/memory.h" + +#include + +#ifdef CFD_ENABLE_OPENMP + +#include + +/* ============================================================================ + * JACOBI CONTEXT + * ============================================================================ */ + +typedef struct { + double dx2; /* dx^2 */ + double dy2; /* dy^2 */ + double inv_dz2; /* 1/dz^2 (0 for 2D) */ + double inv_factor; /* 1 / (2 * (1/dx^2 + 1/dy^2 + inv_dz2)) */ + size_t stride_z; /* nx*ny for 3D, 0 for 2D */ + size_t k_start; /* first interior k index */ + size_t k_end; /* one-past-last interior k index */ + size_t nz; + int initialized; +} jacobi_omp_context_t; + +/* ============================================================================ + * JACOBI OMP IMPLEMENTATION + * ============================================================================ */ + +static cfd_status_t jacobi_omp_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; + + jacobi_omp_context_t* ctx = (jacobi_omp_context_t*)cfd_calloc(1, sizeof(jacobi_omp_context_t)); + if (!ctx) { + return CFD_ERROR_NOMEM; + } + + ctx->dx2 = dx * dx; + ctx->dy2 = dy * dy; + ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); + ctx->nz = nz; + poisson_solver_compute_3d_bounds(nz, nx, ny, + &ctx->stride_z, &ctx->k_start, &ctx->k_end); + + double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2); + ctx->inv_factor = 1.0 / factor; + ctx->initialized = 1; + + solver->context = ctx; + return CFD_SUCCESS; +} + +static void jacobi_omp_destroy(poisson_solver_t* solver) { + if (solver && solver->context) { + cfd_free(solver->context); + solver->context = NULL; + } +} + +static cfd_status_t jacobi_omp_iterate( + poisson_solver_t* solver, + double* x, + double* x_temp, + const double* rhs, + double* residual) +{ + if (!x_temp) { + return CFD_ERROR_INVALID; /* Jacobi requires temp buffer */ + } + + jacobi_omp_context_t* ctx = (jacobi_omp_context_t*)solver->context; + size_t nx = solver->nx; + size_t ny = solver->ny; + double dx2 = ctx->dx2; + double dy2 = ctx->dy2; + double inv_dz2 = ctx->inv_dz2; + double inv_factor = ctx->inv_factor; + size_t stride_z = ctx->stride_z; + + const double* p_old = x; + double* p_new = x_temp; + + /* Jacobi update: reads from p_old, writes to p_new (no in-place dependency) */ + for (size_t k = ctx->k_start; k < ctx->k_end; k++) { + int j; +#pragma omp parallel for schedule(static) + for (j = 1; j < (int)ny - 1; j++) { + for (int i = 1; i < (int)nx - 1; i++) { + size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); + + double p_result = -(rhs[idx] + - (p_old[idx + 1] + p_old[idx - 1]) / dx2 + - (p_old[idx + nx] + p_old[idx - nx]) / dy2 + - (p_old[idx + stride_z] + p_old[idx - stride_z]) * inv_dz2 + ) * inv_factor; + + p_new[idx] = p_result; + } + } + } + + /* Copy result back to x */ + memcpy(x, x_temp, nx * ny * ctx->nz * sizeof(double)); + + /* Apply boundary conditions */ + poisson_solver_apply_bc(solver, x); + + /* Compute residual if requested */ + if (residual) { + *residual = poisson_solver_compute_residual(solver, x, rhs); + } + + return CFD_SUCCESS; +} + +/* ============================================================================ + * FACTORY FUNCTION + * ============================================================================ */ + +poisson_solver_t* create_jacobi_omp_solver(void) { + poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t)); + if (!solver) { + return NULL; + } + + solver->name = POISSON_SOLVER_TYPE_JACOBI_OMP; + solver->description = "Jacobi iteration (OpenMP)"; + solver->method = POISSON_METHOD_JACOBI; + solver->backend = POISSON_BACKEND_OMP; + solver->params = poisson_solver_params_default(); + solver->params.max_iterations = 2000; /* Jacobi needs more iterations */ + solver->params.check_interval = 10; /* Check less frequently for speed */ + + solver->init = jacobi_omp_init; + solver->destroy = jacobi_omp_destroy; + solver->solve = NULL; /* Use common solve loop */ + solver->iterate = jacobi_omp_iterate; + solver->apply_bc = NULL; /* Use default Neumann */ + + return solver; +} + +#endif /* CFD_ENABLE_OPENMP */ diff --git a/tests/math/test_bicgstab.c b/tests/math/test_bicgstab.c index 83e1ffdc..5e6e5bdf 100644 --- a/tests/math/test_bicgstab.c +++ b/tests/math/test_bicgstab.c @@ -539,18 +539,25 @@ void test_bicgstab_dirichlet(void) { * ============================================================================ */ /** - * Test that unsupported backends return NULL + * Test backend availability contract: create() returns a solver exactly when + * the requested backend was compiled in, and NULL otherwise (no silent + * fallbacks). CFD_ENABLE_OPENMP / CFD_HAS_CUDA are library-compile macros not + * visible in this test TU, so query the backend at runtime via + * backend_available(), which tracks the same compile-time support. */ -void test_bicgstab_unsupported_backend(void) { - /* OMP backend not yet implemented for BiCGSTAB */ +void test_bicgstab_backend_availability(void) { + /* OMP backend: available whenever the library was built with OpenMP. */ poisson_solver_t* solver = poisson_solver_create( POISSON_METHOD_BICGSTAB, POISSON_BACKEND_OMP); - TEST_ASSERT_NULL(solver); + if (poisson_solver_backend_available(POISSON_BACKEND_OMP)) { + TEST_ASSERT_NOT_NULL(solver); + poisson_solver_destroy(solver); + } else { + TEST_ASSERT_NULL(solver); + } /* GPU backend: create() returns the factory solver whenever the library was - * built with CUDA (the device check happens in init, not create). CFD_HAS_CUDA - * is a library-compile macro not visible in this test TU, so query the backend - * at runtime — backend_available() tracks the same compile-time CUDA support. */ + * built with CUDA (the device check happens in init, not create). */ solver = poisson_solver_create( POISSON_METHOD_BICGSTAB, POISSON_BACKEND_GPU); if (poisson_solver_backend_available(POISSON_BACKEND_GPU)) { @@ -588,7 +595,7 @@ int main(void) { RUN_TEST(test_bicgstab_dirichlet); /* Error handling tests */ - RUN_TEST(test_bicgstab_unsupported_backend); + RUN_TEST(test_bicgstab_backend_availability); RUN_TEST(test_bicgstab_destroy_null); return UNITY_END(); diff --git a/tests/math/test_omp_consistency.c b/tests/math/test_omp_consistency.c index 4fbacb3d..730c1b59 100644 --- a/tests/math/test_omp_consistency.c +++ b/tests/math/test_omp_consistency.c @@ -398,10 +398,209 @@ void test_gmres_omp_vs_scalar(void) { cfd_free(rhs); } +/** + * Test: Jacobi OMP vs Scalar Consistency + * + * Jacobi's interior update is a pure element-wise map (reads the previous + * iterate, writes the next) with no reduction, and both backends drive + * convergence with the same shared scalar residual check — so the OMP and + * scalar solutions are effectively bit-identical. Jacobi converges slowly, so a + * larger iteration budget is needed at 33x33 (~2900 iterations). + */ +void test_jacobi_omp_vs_scalar(void) { + double dx = (XMAX - XMIN) / (NX - 1); + double dy = (YMAX - YMIN) / (NY - 1); + size_t n = NX * NY; + + double* x_scalar = (double*)cfd_calloc(n, sizeof(double)); + double* x_omp = (double*)cfd_calloc(n, sizeof(double)); + double* x_temp = (double*)cfd_calloc(n, sizeof(double)); + double* rhs = (double*)cfd_calloc(n, sizeof(double)); + TEST_ASSERT_NOT_NULL(x_scalar); + TEST_ASSERT_NOT_NULL(x_omp); + TEST_ASSERT_NOT_NULL(x_temp); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, NX, NY, dx, dy); + + poisson_solver_t* solver_scalar = poisson_solver_create( + POISSON_METHOD_JACOBI, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver_scalar); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = 5000; /* Jacobi needs ~2900 iters at 33x33 */ + + cfd_status_t status = poisson_solver_init(solver_scalar, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_scalar = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_scalar, x_scalar, x_temp, rhs, &stats_scalar); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_scalar.status); + + /* OMP backend requires OpenMP - skip at runtime if unavailable */ + if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); + return; + } + + poisson_solver_t* solver_omp = poisson_solver_create( + POISSON_METHOD_JACOBI, POISSON_BACKEND_OMP); + + if (!solver_omp) { + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(x_temp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("OMP Jacobi solver creation failed; backend may be unavailable"); + return; + } + + status = poisson_solver_init(solver_omp, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_omp = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_omp, x_omp, x_temp, rhs, &stats_omp); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_omp.status); + + double l2_diff = 0.0; + size_t count = 0; + for (size_t j = 1; j < NY - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + size_t idx = IDX_2D(i, j, NX); + double diff = x_scalar[idx] - x_omp[idx]; + l2_diff += diff * diff; + count++; + } + } + l2_diff = sqrt(l2_diff / count); + + /* Jacobi's element-wise update carries no reduction, so OMP and scalar are + * effectively bit-identical. */ + TEST_ASSERT_DOUBLE_WITHIN(1.0e-9, 0.0, l2_diff); + + int iter_diff = abs((int)stats_scalar.iterations - (int)stats_omp.iterations); + TEST_ASSERT_LESS_OR_EQUAL(2, iter_diff); + + poisson_solver_destroy(solver_scalar); + poisson_solver_destroy(solver_omp); + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(x_temp); + cfd_free(rhs); +} + +/** + * Test: BiCGSTAB OMP vs Scalar Consistency + * + * The BiCGSTAB solve loop and every per-element vector update are identical + * across backends; only the dot-product reductions accumulate in a different + * order across threads. Those reductions feed the scalar coefficients + * (alpha/beta/omega), so BiCGSTAB is more path-sensitive than CG — a relaxed + * tolerance (matching the Red-Black SOR consistency test) is used, verifying + * both backends converge to the same solution within solver tolerance. + */ +void test_bicgstab_omp_vs_scalar(void) { + double dx = (XMAX - XMIN) / (NX - 1); + double dy = (YMAX - YMIN) / (NY - 1); + size_t n = NX * NY; + + double* x_scalar = (double*)cfd_calloc(n, sizeof(double)); + double* x_omp = (double*)cfd_calloc(n, sizeof(double)); + double* rhs = (double*)cfd_calloc(n, sizeof(double)); + TEST_ASSERT_NOT_NULL(x_scalar); + TEST_ASSERT_NOT_NULL(x_omp); + TEST_ASSERT_NOT_NULL(rhs); + + init_sinusoidal_rhs(rhs, NX, NY, dx, dy); + + poisson_solver_t* solver_scalar = poisson_solver_create( + POISSON_METHOD_BICGSTAB, POISSON_BACKEND_SCALAR); + TEST_ASSERT_NOT_NULL(solver_scalar); + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = TOLERANCE; + params.max_iterations = 1000; + + cfd_status_t status = poisson_solver_init(solver_scalar, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_scalar = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_scalar, x_scalar, NULL, rhs, &stats_scalar); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_scalar.status); + + /* OMP backend requires OpenMP - skip at runtime if unavailable */ + if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); + return; + } + + poisson_solver_t* solver_omp = poisson_solver_create( + POISSON_METHOD_BICGSTAB, POISSON_BACKEND_OMP); + + if (!solver_omp) { + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(rhs); + poisson_solver_destroy(solver_scalar); + TEST_IGNORE_MESSAGE("OMP BiCGSTAB solver creation failed; backend may be unavailable"); + return; + } + + status = poisson_solver_init(solver_omp, NX, NY, 1, dx, dy, 0.0, ¶ms); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + + poisson_solver_stats_t stats_omp = poisson_solver_stats_default(); + status = poisson_solver_solve(solver_omp, x_omp, NULL, rhs, &stats_omp); + TEST_ASSERT_EQUAL(CFD_SUCCESS, status); + TEST_ASSERT_EQUAL(POISSON_CONVERGED, stats_omp.status); + + double l2_diff = 0.0; + size_t count = 0; + for (size_t j = 1; j < NY - 1; j++) { + for (size_t i = 1; i < NX - 1; i++) { + size_t idx = IDX_2D(i, j, NX); + double diff = x_scalar[idx] - x_omp[idx]; + l2_diff += diff * diff; + count++; + } + } + l2_diff = sqrt(l2_diff / count); + + /* Both backends converge to the same solution within solver tolerance; + * reduction-order differences in alpha/beta/omega make BiCGSTAB more + * path-sensitive than CG, so the tolerance is relaxed (as for RB-SOR). */ + TEST_ASSERT_DOUBLE_WITHIN(1.0e-6, 0.0, l2_diff); + + int iter_diff = abs((int)stats_scalar.iterations - (int)stats_omp.iterations); + TEST_ASSERT_LESS_OR_EQUAL(5, iter_diff); + + poisson_solver_destroy(solver_scalar); + poisson_solver_destroy(solver_omp); + cfd_free(x_scalar); + cfd_free(x_omp); + cfd_free(rhs); +} + int main(void) { UNITY_BEGIN(); RUN_TEST(test_cg_omp_vs_scalar); RUN_TEST(test_redblack_omp_vs_scalar); RUN_TEST(test_gmres_omp_vs_scalar); + RUN_TEST(test_jacobi_omp_vs_scalar); + RUN_TEST(test_bicgstab_omp_vs_scalar); return UNITY_END(); } From 3bad5d3bfdf98a86d3b5fe2a0e3fd6c57c243dbf Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 12 Sep 2026 11:23:43 +0300 Subject: [PATCH 2/6] Reject grids above INT_MAX in Jacobi OMP init The interior sweep casts nx/ny to int for its OpenMP loop bounds; a larger dimension overflows the cast and the sweep silently does no work. Fail init with CFD_ERROR_LIMIT_EXCEEDED so the solver never runs on unsupported dims. --- .../linear/omp/linear_solver_jacobi_omp.c | 6 ++++ tests/math/test_omp_consistency.c | 30 +++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c b/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c index 6826ed84..c14b0a55 100644 --- a/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c @@ -47,6 +47,12 @@ static cfd_status_t jacobi_omp_init( { (void)params; + /* The interior sweep uses int OpenMP loop bounds; larger dims would overflow the cast */ + if (nx > (size_t)INT_MAX || ny > (size_t)INT_MAX) { + cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop"); + return CFD_ERROR_LIMIT_EXCEEDED; + } + jacobi_omp_context_t* ctx = (jacobi_omp_context_t*)cfd_calloc(1, sizeof(jacobi_omp_context_t)); if (!ctx) { return CFD_ERROR_NOMEM; diff --git a/tests/math/test_omp_consistency.c b/tests/math/test_omp_consistency.c index 730c1b59..3c85e219 100644 --- a/tests/math/test_omp_consistency.c +++ b/tests/math/test_omp_consistency.c @@ -15,6 +15,7 @@ #include "cfd/solvers/poisson_solver.h" #include "cfd/core/memory.h" #include "cfd/core/indexing.h" +#include #include #include @@ -595,6 +596,34 @@ void test_bicgstab_omp_vs_scalar(void) { cfd_free(rhs); } +/** + * Test: Jacobi OMP init rejects dimensions above INT_MAX + * + * The interior sweep uses int OpenMP loop bounds, so init must refuse a grid + * dimension that would overflow the cast instead of running a truncated sweep. + */ +void test_jacobi_omp_init_rejects_dims_above_int_max(void) { + if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { + TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); + return; + } + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_JACOBI, POISSON_BACKEND_OMP); + if (!solver) { + TEST_IGNORE_MESSAGE("OMP Jacobi solver creation failed; backend may be unavailable"); + return; + } + + size_t too_large = (size_t)INT_MAX + 1; + TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, + poisson_solver_init(solver, too_large, NY, 1, 1.0, 1.0, 0.0, NULL)); + TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, + poisson_solver_init(solver, NX, too_large, 1, 1.0, 1.0, 0.0, NULL)); + + poisson_solver_destroy(solver); +} + int main(void) { UNITY_BEGIN(); RUN_TEST(test_cg_omp_vs_scalar); @@ -602,5 +631,6 @@ int main(void) { RUN_TEST(test_gmres_omp_vs_scalar); RUN_TEST(test_jacobi_omp_vs_scalar); RUN_TEST(test_bicgstab_omp_vs_scalar); + RUN_TEST(test_jacobi_omp_init_rejects_dims_above_int_max); return UNITY_END(); } From 2d7fa742471de6fa541796f31131de9563941bab Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 12 Sep 2026 11:25:03 +0300 Subject: [PATCH 3/6] Reject grids above INT_MAX in BiCGSTAB OMP init bicgstab_size_to_int() maps an nx/ny above INT_MAX to 0, so every OMP primitive skips its sweep, the initial residual reads zero, and the solver reports false convergence. Fail init with CFD_ERROR_LIMIT_EXCEEDED before allocating the work vectors. --- .../linear/omp/linear_solver_bicgstab_omp.c | 7 +++++ tests/math/test_omp_consistency.c | 29 +++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c index b4657638..2430a6ba 100644 --- a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c @@ -286,6 +286,13 @@ static cfd_status_t bicgstab_omp_init( { (void)params; + /* The primitives use int OpenMP loop bounds; bicgstab_size_to_int() maps larger + * dims to 0, which would turn every sweep into a no-op and fake convergence */ + if (nx > (size_t)INT_MAX || ny > (size_t)INT_MAX) { + cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop"); + return CFD_ERROR_LIMIT_EXCEEDED; + } + bicgstab_omp_context_t* ctx = (bicgstab_omp_context_t*)cfd_calloc(1, sizeof(bicgstab_omp_context_t)); if (!ctx) { return CFD_ERROR_NOMEM; diff --git a/tests/math/test_omp_consistency.c b/tests/math/test_omp_consistency.c index 3c85e219..89e82537 100644 --- a/tests/math/test_omp_consistency.c +++ b/tests/math/test_omp_consistency.c @@ -624,6 +624,34 @@ void test_jacobi_omp_init_rejects_dims_above_int_max(void) { poisson_solver_destroy(solver); } +/** + * Test: BiCGSTAB OMP init rejects dimensions above INT_MAX + * + * The OMP primitives use int loop bounds and map an overflowing dimension to 0, + * which would skip every sweep and report false convergence — init must refuse. + */ +void test_bicgstab_omp_init_rejects_dims_above_int_max(void) { + if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { + TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); + return; + } + + poisson_solver_t* solver = poisson_solver_create( + POISSON_METHOD_BICGSTAB, POISSON_BACKEND_OMP); + if (!solver) { + TEST_IGNORE_MESSAGE("OMP BiCGSTAB solver creation failed; backend may be unavailable"); + return; + } + + size_t too_large = (size_t)INT_MAX + 1; + TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, + poisson_solver_init(solver, too_large, NY, 1, 1.0, 1.0, 0.0, NULL)); + TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, + poisson_solver_init(solver, NX, too_large, 1, 1.0, 1.0, 0.0, NULL)); + + poisson_solver_destroy(solver); +} + int main(void) { UNITY_BEGIN(); RUN_TEST(test_cg_omp_vs_scalar); @@ -632,5 +660,6 @@ int main(void) { RUN_TEST(test_jacobi_omp_vs_scalar); RUN_TEST(test_bicgstab_omp_vs_scalar); RUN_TEST(test_jacobi_omp_init_rejects_dims_above_int_max); + RUN_TEST(test_bicgstab_omp_init_rejects_dims_above_int_max); return UNITY_END(); } From 1e9a9df9dff9e7bb1a801dff3c55a58cd99e16c3 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 12 Sep 2026 11:35:37 +0300 Subject: [PATCH 4/6] Use renamed poisson_solver_size_to_int in BiCGSTAB OMP Master's #206 renamed bicgstab_size_to_int, so once master is merged the OMP BiCGSTAB primitives call an undeclared function and every CI build fails to link. --- .../linear/omp/linear_solver_bicgstab_omp.c | 42 +++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c index 2430a6ba..2fbe55ff 100644 --- a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c @@ -57,8 +57,8 @@ static double dot_product_omp(const double* a, const double* b, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { double sum = 0.0; - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -77,8 +77,8 @@ static double dot_product_omp(const double* a, const double* b, static void axpy_omp(double alpha, const double* x, double* y, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -103,8 +103,8 @@ static void apply_laplacian_omp(const double* p, double* Ap, size_t k_start, size_t k_end, size_t stride_z) { double dx2_inv = 1.0 / dx2; double dy2_inv = 1.0 / dy2; - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -132,8 +132,8 @@ static void compute_residual_omp(const double* x, const double* rhs, double* r, size_t k_start, size_t k_end, size_t stride_z) { double dx2_inv = 1.0 / dx2; double dy2_inv = 1.0 / dy2; - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -155,8 +155,8 @@ static void compute_residual_omp(const double* x, const double* rhs, double* r, static void copy_vector_omp(const double* src, double* dst, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -173,8 +173,8 @@ static void copy_vector_omp(const double* src, double* dst, /* v = 0 (interior points only) */ static void zero_vector_omp(double* v, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -199,8 +199,8 @@ static void update_p_omp(double* p, const double* r, const double* v, double beta, double omega, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -219,8 +219,8 @@ static void update_s_omp(double* s, const double* r, const double* v, double alpha, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -239,8 +239,8 @@ static void update_x_omp(double* x, const double* p, const double* s, double alpha, double omega, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -259,8 +259,8 @@ static void update_r_omp(double* r, const double* s, const double* t, double omega, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = bicgstab_size_to_int(ny); - int nx_int = bicgstab_size_to_int(nx); + int ny_int = poisson_solver_size_to_int(ny); + int nx_int = poisson_solver_size_to_int(nx); for (size_t k = k_start; k < k_end; k++) { int j; @@ -286,7 +286,7 @@ static cfd_status_t bicgstab_omp_init( { (void)params; - /* The primitives use int OpenMP loop bounds; bicgstab_size_to_int() maps larger + /* The primitives use int OpenMP loop bounds; poisson_solver_size_to_int() maps larger * dims to 0, which would turn every sweep into a no-op and fake convergence */ if (nx > (size_t)INT_MAX || ny > (size_t)INT_MAX) { cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop"); From ee868cb4730e320c9201785bffc962db42b9bf35 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 12 Sep 2026 11:37:35 +0300 Subject: [PATCH 5/6] Reuse shared OMP primitives in BiCGSTAB OMP Master's linear_solver_primitives_omp.h already provides dot_product, axpy, apply_laplacian, compute_residual and copy_vector for the OMP Krylov solvers; private copies here would let the implementations drift apart. --- .../linear/omp/linear_solver_bicgstab_omp.c | 125 +----------------- 1 file changed, 7 insertions(+), 118 deletions(-) diff --git a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c index 2fbe55ff..0d803943 100644 --- a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c @@ -25,6 +25,8 @@ #include +#include "linear_solver_primitives_omp.h" + /* ============================================================================ * BICGSTAB CONTEXT * ============================================================================ */ @@ -51,125 +53,12 @@ typedef struct { /* ============================================================================ * OMP-PARALLELIZED PRIMITIVES + * + * dot_product_omp, axpy_omp, apply_laplacian_omp, compute_residual_omp and + * copy_vector_omp come from linear_solver_primitives_omp.h (shared with CG and + * GMRES OMP). * ============================================================================ */ -static double dot_product_omp(const double* a, const double* b, - size_t nx, size_t ny, - size_t k_start, size_t k_end, size_t stride_z) { - double sum = 0.0; - int ny_int = poisson_solver_size_to_int(ny); - int nx_int = poisson_solver_size_to_int(nx); - - for (size_t k = k_start; k < k_end; k++) { - int j; -#pragma omp parallel for schedule(static) reduction(+:sum) - for (j = 1; j < ny_int - 1; j++) { - for (int i = 1; i < nx_int - 1; i++) { - size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); - sum += a[idx] * b[idx]; - } - } - } - return sum; -} - -/* y = y + alpha * x (interior points only) */ -static void axpy_omp(double alpha, const double* x, double* y, - size_t nx, size_t ny, - size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = poisson_solver_size_to_int(ny); - int nx_int = poisson_solver_size_to_int(nx); - - for (size_t k = k_start; k < k_end; k++) { - int j; -#pragma omp parallel for schedule(static) - for (j = 1; j < ny_int - 1; j++) { - for (int i = 1; i < nx_int - 1; i++) { - size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); - y[idx] += alpha * x[idx]; - } - } - } -} - -/* - * Apply negative Laplacian operator: Ap = -nabla^2(p) - * For Poisson equation nabla^2(x) = rhs we solve -nabla^2(x) = -rhs, which is - * SPD with positive eigenvalues. - */ -static void apply_laplacian_omp(const double* p, double* Ap, - size_t nx, size_t ny, - double dx2, double dy2, double inv_dz2, - size_t k_start, size_t k_end, size_t stride_z) { - double dx2_inv = 1.0 / dx2; - double dy2_inv = 1.0 / dy2; - int ny_int = poisson_solver_size_to_int(ny); - int nx_int = poisson_solver_size_to_int(nx); - - for (size_t k = k_start; k < k_end; k++) { - int j; -#pragma omp parallel for schedule(static) - for (j = 1; j < ny_int - 1; j++) { - for (int i = 1; i < nx_int - 1; i++) { - size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); - double laplacian = - (p[idx + 1] - 2.0 * p[idx] + p[idx - 1]) * dx2_inv - + (p[idx + nx] - 2.0 * p[idx] + p[idx - nx]) * dy2_inv - + (p[idx + stride_z] + p[idx - stride_z] - 2.0 * p[idx]) * inv_dz2; - Ap[idx] = -laplacian; - } - } - } -} - -/* - * Compute initial residual: r = b - A*x - * For our formulation: r = -rhs - (-nabla^2 x) = -rhs + nabla^2 x - */ -static void compute_residual_omp(const double* x, const double* rhs, double* r, - size_t nx, size_t ny, - double dx2, double dy2, double inv_dz2, - size_t k_start, size_t k_end, size_t stride_z) { - double dx2_inv = 1.0 / dx2; - double dy2_inv = 1.0 / dy2; - int ny_int = poisson_solver_size_to_int(ny); - int nx_int = poisson_solver_size_to_int(nx); - - for (size_t k = k_start; k < k_end; k++) { - int j; -#pragma omp parallel for schedule(static) - for (j = 1; j < ny_int - 1; j++) { - for (int i = 1; i < nx_int - 1; i++) { - size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); - double laplacian = - (x[idx + 1] - 2.0 * x[idx] + x[idx - 1]) * dx2_inv - + (x[idx + nx] - 2.0 * x[idx] + x[idx - nx]) * dy2_inv - + (x[idx + stride_z] + x[idx - stride_z] - 2.0 * x[idx]) * inv_dz2; - r[idx] = -rhs[idx] + laplacian; - } - } - } -} - -/* dst = src (interior points only) */ -static void copy_vector_omp(const double* src, double* dst, - size_t nx, size_t ny, - size_t k_start, size_t k_end, size_t stride_z) { - int ny_int = poisson_solver_size_to_int(ny); - int nx_int = poisson_solver_size_to_int(nx); - - for (size_t k = k_start; k < k_end; k++) { - int j; -#pragma omp parallel for schedule(static) - for (j = 1; j < ny_int - 1; j++) { - for (int i = 1; i < nx_int - 1; i++) { - size_t idx = k * stride_z + IDX_2D((size_t)i, (size_t)j, nx); - dst[idx] = src[idx]; - } - } - } -} - /* v = 0 (interior points only) */ static void zero_vector_omp(double* v, size_t nx, size_t ny, size_t k_start, size_t k_end, size_t stride_z) { @@ -191,7 +80,7 @@ static void zero_vector_omp(double* v, size_t nx, size_t ny, /* * Fused per-iteration vector updates. Each preserves the scalar reference's * exact per-element floating-point operation order (only the dot-product - * reductions above differ across threads), so OMP and scalar agree tightly. + * reductions differ across threads), so OMP and scalar agree tightly. */ /* p = r + beta * (p - omega * v) */ From 3d24962ff5ebe345c2ae68faeab3327961e367bc Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 12 Sep 2026 11:40:46 +0300 Subject: [PATCH 6/6] Validate OMP Jacobi/BiCGSTAB grid size with the shared helper Master's poisson_solver_validate_grid_size also rejects grids whose nx*ny*nz wraps size_t, which the hand-rolled INT_MAX checks missed: BiCGSTAB OMP would allocate undersized work vectors and Jacobi OMP would memcpy a wrapped size. Cover both in the shared oversized-grid table test instead of one-off tests. --- CHANGELOG.md | 6 +- .../linear/omp/linear_solver_bicgstab_omp.c | 12 ++-- .../linear/omp/linear_solver_jacobi_omp.c | 16 ++--- tests/math/test_omp_consistency.c | 59 ------------------- tests/solvers/test_linear_solver.c | 2 + 5 files changed, 21 insertions(+), 74 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 400d729f..30aaf82f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,11 +13,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 tier (previously CG, GMRES, and Red-Black SOR). Both are selectable via `poisson_solver_create(method, POISSON_BACKEND_OMP)`; per-element updates are identical to the scalar reference, so results match within reduction rounding - (verified by scalar-vs-OMP consistency tests). Plain lexicographic SOR remains + (verified by scalar-vs-OMP consistency tests). Both reject grids whose `nx` or + `ny` exceeds `INT_MAX`, or whose `nx*ny*nz` overflows `size_t`, with + `CFD_ERROR_LIMIT_EXCEEDED` at init. Plain lexicographic SOR remains scalar/SIMD-only — its parallel form is the existing Red-Black SOR OMP solver (`lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c`, `lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c`, - `tests/math/test_omp_consistency.c`). + `tests/math/test_omp_consistency.c`, `tests/solvers/test_linear_solver.c`). - **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 diff --git a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c index 0d803943..0902e939 100644 --- a/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_bicgstab_omp.c @@ -175,11 +175,12 @@ static cfd_status_t bicgstab_omp_init( { (void)params; - /* The primitives use int OpenMP loop bounds; poisson_solver_size_to_int() maps larger - * dims to 0, which would turn every sweep into a no-op and fake convergence */ - if (nx > (size_t)INT_MAX || ny > (size_t)INT_MAX) { - cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop"); - return CFD_ERROR_LIMIT_EXCEEDED; + /* The shared primitives loop over int bounds, which collapse to empty loops + * for oversized dimensions; reject those grids before allocating. */ + size_t n = 0; + cfd_status_t size_status = poisson_solver_validate_grid_size(nx, ny, nz, 1, &n); + if (size_status != CFD_SUCCESS) { + return size_status; } bicgstab_omp_context_t* ctx = (bicgstab_omp_context_t*)cfd_calloc(1, sizeof(bicgstab_omp_context_t)); @@ -193,7 +194,6 @@ static cfd_status_t bicgstab_omp_init( poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end); /* Allocate working vectors */ - size_t n = nx * ny * nz; ctx->r = (double*)cfd_calloc(n, sizeof(double)); ctx->r_hat = (double*)cfd_calloc(n, sizeof(double)); ctx->p = (double*)cfd_calloc(n, sizeof(double)); diff --git a/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c b/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c index c14b0a55..6ad78c21 100644 --- a/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c +++ b/lib/src/solvers/linear/omp/linear_solver_jacobi_omp.c @@ -31,7 +31,7 @@ typedef struct { size_t stride_z; /* nx*ny for 3D, 0 for 2D */ size_t k_start; /* first interior k index */ size_t k_end; /* one-past-last interior k index */ - size_t nz; + size_t n; /* nx*ny*nz, validated at init */ int initialized; } jacobi_omp_context_t; @@ -47,10 +47,12 @@ static cfd_status_t jacobi_omp_init( { (void)params; - /* The interior sweep uses int OpenMP loop bounds; larger dims would overflow the cast */ - if (nx > (size_t)INT_MAX || ny > (size_t)INT_MAX) { - cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop"); - return CFD_ERROR_LIMIT_EXCEEDED; + /* The sweep loops over int bounds and copies nx*ny*nz doubles; oversized + * grids would empty the loops or wrap the copy size, so reject them. */ + size_t n = 0; + cfd_status_t size_status = poisson_solver_validate_grid_size(nx, ny, nz, 1, &n); + if (size_status != CFD_SUCCESS) { + return size_status; } jacobi_omp_context_t* ctx = (jacobi_omp_context_t*)cfd_calloc(1, sizeof(jacobi_omp_context_t)); @@ -61,7 +63,7 @@ static cfd_status_t jacobi_omp_init( ctx->dx2 = dx * dx; ctx->dy2 = dy * dy; ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); - ctx->nz = nz; + ctx->n = n; poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end); @@ -123,7 +125,7 @@ static cfd_status_t jacobi_omp_iterate( } /* Copy result back to x */ - memcpy(x, x_temp, nx * ny * ctx->nz * sizeof(double)); + memcpy(x, x_temp, ctx->n * sizeof(double)); /* Apply boundary conditions */ poisson_solver_apply_bc(solver, x); diff --git a/tests/math/test_omp_consistency.c b/tests/math/test_omp_consistency.c index 89e82537..730c1b59 100644 --- a/tests/math/test_omp_consistency.c +++ b/tests/math/test_omp_consistency.c @@ -15,7 +15,6 @@ #include "cfd/solvers/poisson_solver.h" #include "cfd/core/memory.h" #include "cfd/core/indexing.h" -#include #include #include @@ -596,62 +595,6 @@ void test_bicgstab_omp_vs_scalar(void) { cfd_free(rhs); } -/** - * Test: Jacobi OMP init rejects dimensions above INT_MAX - * - * The interior sweep uses int OpenMP loop bounds, so init must refuse a grid - * dimension that would overflow the cast instead of running a truncated sweep. - */ -void test_jacobi_omp_init_rejects_dims_above_int_max(void) { - if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { - TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); - return; - } - - poisson_solver_t* solver = poisson_solver_create( - POISSON_METHOD_JACOBI, POISSON_BACKEND_OMP); - if (!solver) { - TEST_IGNORE_MESSAGE("OMP Jacobi solver creation failed; backend may be unavailable"); - return; - } - - size_t too_large = (size_t)INT_MAX + 1; - TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, - poisson_solver_init(solver, too_large, NY, 1, 1.0, 1.0, 0.0, NULL)); - TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, - poisson_solver_init(solver, NX, too_large, 1, 1.0, 1.0, 0.0, NULL)); - - poisson_solver_destroy(solver); -} - -/** - * Test: BiCGSTAB OMP init rejects dimensions above INT_MAX - * - * The OMP primitives use int loop bounds and map an overflowing dimension to 0, - * which would skip every sweep and report false convergence — init must refuse. - */ -void test_bicgstab_omp_init_rejects_dims_above_int_max(void) { - if (!poisson_solver_backend_available(POISSON_BACKEND_OMP)) { - TEST_IGNORE_MESSAGE("OMP backend not available on this platform"); - return; - } - - poisson_solver_t* solver = poisson_solver_create( - POISSON_METHOD_BICGSTAB, POISSON_BACKEND_OMP); - if (!solver) { - TEST_IGNORE_MESSAGE("OMP BiCGSTAB solver creation failed; backend may be unavailable"); - return; - } - - size_t too_large = (size_t)INT_MAX + 1; - TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, - poisson_solver_init(solver, too_large, NY, 1, 1.0, 1.0, 0.0, NULL)); - TEST_ASSERT_EQUAL(CFD_ERROR_LIMIT_EXCEEDED, - poisson_solver_init(solver, NX, too_large, 1, 1.0, 1.0, 0.0, NULL)); - - poisson_solver_destroy(solver); -} - int main(void) { UNITY_BEGIN(); RUN_TEST(test_cg_omp_vs_scalar); @@ -659,7 +602,5 @@ int main(void) { RUN_TEST(test_gmres_omp_vs_scalar); RUN_TEST(test_jacobi_omp_vs_scalar); RUN_TEST(test_bicgstab_omp_vs_scalar); - RUN_TEST(test_jacobi_omp_init_rejects_dims_above_int_max); - RUN_TEST(test_bicgstab_omp_init_rejects_dims_above_int_max); return UNITY_END(); } diff --git a/tests/solvers/test_linear_solver.c b/tests/solvers/test_linear_solver.c index 8c320a69..20f82a55 100644 --- a/tests/solvers/test_linear_solver.c +++ b/tests/solvers/test_linear_solver.c @@ -1104,6 +1104,8 @@ void test_poisson_init_rejects_oversized_grid(void) { poisson_solver_backend_t backend; } solvers[] = { { POISSON_METHOD_CG, POISSON_BACKEND_OMP }, + { POISSON_METHOD_JACOBI, POISSON_BACKEND_OMP }, + { POISSON_METHOD_BICGSTAB, POISSON_BACKEND_OMP }, { POISSON_METHOD_BICGSTAB, POISSON_BACKEND_SIMD }, }; const size_t dims[][3] = {