Add Jacobi and BiCGSTAB OpenMP linear solvers

CHANGELOG.md

@@ -9,6 +9,17 @@ 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). 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/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

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

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
 )
 

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"

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();

lib/src/solvers/linear/linear_solver_internal.h

@@ -42,8 +42,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 */