Add GPU plain SOR (Block SOR) linear solver, closing last backend gap

CMakeLists.txt

@@ -319,6 +319,7 @@ if(BUILD_TESTS)
     add_executable(test_finite_differences_3d tests/math/test_finite_differences_3d.c)
     add_executable(test_poisson_accuracy tests/math/test_poisson_accuracy.c)
     add_executable(test_poisson_jacobi_gpu tests/math/test_poisson_jacobi_gpu.c)
+    add_executable(test_poisson_sor_gpu tests/math/test_poisson_sor_gpu.c)
     add_executable(test_poisson_3d tests/math/test_poisson_3d.c)
     add_executable(test_laplacian_accuracy tests/math/test_laplacian_accuracy.c)
     add_executable(test_linear_solver_convergence tests/math/test_linear_solver_convergence.c)
@@ -431,6 +432,7 @@ if(BUILD_TESTS)
     target_include_directories(test_finite_differences_3d PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lib/include)
     target_link_libraries(test_poisson_accuracy PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_poisson_jacobi_gpu PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
+    target_link_libraries(test_poisson_sor_gpu PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_poisson_3d PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_laplacian_accuracy PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_link_libraries(test_linear_solver_convergence PRIVATE CFD::Library unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
@@ -577,6 +579,7 @@ if(BUILD_TESTS)
     add_test(NAME FiniteDifferences3DTest COMMAND test_finite_differences_3d)
     add_test(NAME PoissonAccuracyTest COMMAND test_poisson_accuracy)
     add_test(NAME PoissonJacobiGpuTest COMMAND test_poisson_jacobi_gpu)
+    add_test(NAME PoissonSorGpuTest COMMAND test_poisson_sor_gpu)
     add_test(NAME Poisson3DTest COMMAND test_poisson_3d)
     add_test(NAME LaplacianAccuracyTest COMMAND test_laplacian_accuracy)
     add_test(NAME LinearSolverConvergenceTest COMMAND test_linear_solver_convergence)

ROADMAP.md

@@ -48,7 +48,7 @@ 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 |
 | **Linear Solvers**  | Jacobi         | done | done     | done     | —        | done |
-|                     | SOR            | 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 |
@@ -104,8 +104,9 @@ No-slip, Inlet, Outlet, Symmetry, Moving wall, Time-varying). See CHANGELOG.
 ### 1.2 Linear Solvers (P0)
 
 Implemented: Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB (all with SIMD backends; CG is the
-default Poisson solver for projection methods). GPU standalone Jacobi, CG, Red-Black SOR, and
-BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses on-device CG.
+default Poisson solver for projection methods). GPU standalone Jacobi, CG, Red-Black SOR,
+plain SOR, and BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses
+on-device CG.
 
 **Still needed:**
 
@@ -115,7 +116,7 @@ BiCGSTAB are done and validated vs CPU; `solve_projection_method_gpu` uses on-de
 - [ ] ILU preconditioner
 - [ ] Geometric multigrid
 - [ ] Algebraic multigrid (AMG) — solver and preconditioner (for CG/GMRES/BiCGSTAB)
-- [ ] GPU plain SOR (the one remaining backend gap in the matrix above)
+- [x] GPU plain SOR (Block SOR: per-thread tile sweep, red-black tile coloring, in-place; closes the matrix)
 
 ### 1.3 Numerical Schemes (P1)
 

docs/reference/solvers.md

@@ -125,6 +125,7 @@ p_ij^(k+1) = (1-ω)p_ij^k + (ω/4)(p_i-1,j + p_i+1,j + p_i,j-1 + p_i,j+1 - h²f_
 - Sequential row updates (row j depends on j-1)
 - Optimal ω depends on problem
 - SIMD variant uses Block SOR: processes SIMD_WIDTH consecutive cells per block, with intra-block left-neighbor approximation (see [Block SOR technical note](../technical-notes/block-sor-simd.md))
+- GPU variant also uses Block SOR: each thread sweeps an 8×8 tile sequentially (Gauss-Seidel inside the tile), with red-black *tile* coloring (red pass then black pass, two launches per iteration) so a tile's halo is never written by another tile in the same pass — the update is in-place, race-free, and provably convergent for 0<ω<2
 
 **Convergence Rate:** ρ ≈ 1 - 2πh (with optimal ω)
 
@@ -133,6 +134,7 @@ p_ij^(k+1) = (1-ω)p_ij^k + (ω/4)(p_i-1,j + p_i+1,j + p_i,j-1 + p_i,j+1 - h²f_
 |--------|---------|-------------|
 | `sor_scalar` | Scalar | Sequential Gauss-Seidel + SOR relaxation |
 | `sor_simd` | SIMD | Block SOR (auto-detects AVX2/NEON) |
+| `sor_gpu` | GPU | Block SOR (CUDA; per-thread tile sweep, red-black tile coloring, in-place) |
 
 **Usage:**
 ```c

lib/CMakeLists.txt

@@ -194,6 +194,7 @@ set(CFD_CUDA_SOURCES
     src/solvers/linear/gpu/poisson_solver_cg_gpu.cu
     src/solvers/linear/gpu/poisson_solver_bicgstab_gpu.cu
     src/solvers/linear/gpu/poisson_solver_redblack_sor_gpu.cu
+    src/solvers/linear/gpu/poisson_solver_sor_gpu.cu
 )
 
 #=============================================================================

lib/src/solvers/linear/gpu/poisson_gpu_primitives.cuh

@@ -89,6 +89,61 @@ static __global__ void lin_gpu_kernel_redblack_sweep(double* __restrict__ x,
     }
 }
 
+/**
+ * Block (tiled) SOR sweep: each thread does a sequential Gauss-Seidel/SOR sweep
+ * over one tile_w x tile_h tile of the interior, in row-major order. Tiles are
+ * parallelized by a red-black *tile* coloring (color = (tile_col+tile_row+k)&1):
+ * a single launch updates only tiles of the requested `color`, and every neighbor
+ * tile has the opposite color (adjacency flips parity), so a tile's halo cells are
+ * never written by another thread in the same pass. Running the two colors as two
+ * launches (the launch boundary is the color sync) gives a consistent ordering, so
+ * unlike a Jacobi-coupled block scheme this is provably convergent for 0<omega<2 —
+ * the auto optimal-omega stays stable.
+ *
+ * Updates are in-place: within the tile the left/down neighbors are already swept
+ * (fresh) and the right/up neighbors are not yet swept (still this-iteration's
+ * start values) — exactly the lexicographic SOR ordering. Cross-tile and
+ * z-neighbors belong to opposite-color tiles, stable this pass. At 1x1 tiles this
+ * degenerates to the cell-level Red-Black SOR. Per-cell update matches the SOR
+ * reference:  p_new = (sum_neighbors - rhs)*inv_factor;  x += omega*(p_new - x).
+ */
+static __global__ void lin_gpu_kernel_block_sor_tile_sweep(double* __restrict__ x,
+                                                           const double* __restrict__ rhs,
+                                                           int color, double omega,
+                                                           int tile_w, int tile_h,
+                                                           size_t nx, size_t ny,
+                                                           size_t stride_z, int k_start, int k_end,
+                                                           double inv_dx2, double inv_dy2, double inv_dz2,
+                                                           double inv_factor) {
+    int tile_col = blockIdx.x * blockDim.x + threadIdx.x;
+    int tile_row = blockIdx.y * blockDim.y + threadIdx.y;
+    int i0 = 1 + tile_col * tile_w;  /* tile origin in interior index space [1, n-1) */
+    int j0 = 1 + tile_row * tile_h;
+    if (i0 >= (int)nx - 1 || j0 >= (int)ny - 1)
+        return;
+
+    for (int k = k_start; k <= k_end; k++) {
+        if (((tile_col + tile_row + k) & 1) != color)
+            continue;  /* this tile belongs to the other color on this k-plane */
+        for (int dj = 0; dj < tile_h; dj++) {
+            int j = j0 + dj;
+            if (j >= (int)ny - 1)
+                break;
+            for (int di = 0; di < tile_w; di++) {
+                int i = i0 + di;
+                if (i >= (int)nx - 1)
+                    break;
+                size_t idx = (size_t)k * stride_z + IDX_2D(i, j, nx);
+                double sum = (x[idx - 1] + x[idx + 1]) * inv_dx2
+                           + (x[idx - nx] + x[idx + nx]) * inv_dy2
+                           + (x[idx - stride_z] + x[idx + stride_z]) * inv_dz2;
+                double p_new = (sum - rhs[idx]) * inv_factor;
+                x[idx] += omega * (p_new - x[idx]);
+            }
+        }
+    }
+}
+
 /**
  * Apply the Laplacian operator: out = A*x where A is the 5/7-point Laplacian.
  *   out = sum_neighbors - factor*x