diff --git a/CMakeLists.txt b/CMakeLists.txt index c70c56b6..036bd5cc 100644 --- a/CMakeLists.txt +++ b/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 $<$>:m>) target_link_libraries(test_poisson_jacobi_gpu PRIVATE CFD::Library unity $<$>:m>) + target_link_libraries(test_poisson_sor_gpu PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_poisson_3d PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_laplacian_accuracy PRIVATE CFD::Library unity $<$>:m>) target_link_libraries(test_linear_solver_convergence PRIVATE CFD::Library unity $<$>: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) diff --git a/ROADMAP.md b/ROADMAP.md index 0e2abeef..7dbc1b70 100644 --- a/ROADMAP.md +++ b/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) diff --git a/docs/reference/solvers.md b/docs/reference/solvers.md index f523e0de..0cb6a0c2 100644 --- a/docs/reference/solvers.md +++ b/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 diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index 153f05bd..8300d605 100644 --- a/lib/CMakeLists.txt +++ b/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 ) #============================================================================= diff --git a/lib/src/solvers/linear/gpu/poisson_gpu_primitives.cuh b/lib/src/solvers/linear/gpu/poisson_gpu_primitives.cuh index 9be45e7a..97d52ad7 100644 --- a/lib/src/solvers/linear/gpu/poisson_gpu_primitives.cuh +++ b/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= (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 diff --git a/lib/src/solvers/linear/gpu/poisson_solver_sor_gpu.cu b/lib/src/solvers/linear/gpu/poisson_solver_sor_gpu.cu new file mode 100644 index 00000000..9b8eff81 --- /dev/null +++ b/lib/src/solvers/linear/gpu/poisson_solver_sor_gpu.cu @@ -0,0 +1,305 @@ +/** + * @file poisson_solver_sor_gpu.cu + * @brief Plain SOR Poisson solver — CUDA GPU backend (Block SOR) + * + * Implements the poisson_solver_t interface on the GPU. Like the Jacobi and + * Red-Black GPU backends (and unlike the CPU/SIMD plain-SOR solvers, which expose + * a per-iteration `iterate` driven by the common host solve loop), this backend + * implements `solve` directly: it uploads the RHS and initial guess once, runs the + * full iteration on-device with a relative-residual convergence check, then + * downloads the result. Per-iteration host/device transfers would dominate. + * + * Plain lexicographic SOR is sequential (each cell depends on its already-updated + * left/down neighbors), which the GPU cannot reproduce cheaply. This backend uses + * the same "Block SOR" idea as the AVX2/NEON plain-SOR solvers — each thread does a + * sequential Gauss-Seidel/SOR sweep over a small tile — but parallelizes the tiles + * with a red-black *tile* coloring: each iteration launches the red tile pass then + * the black tile pass (the launch boundary is the color sync). Because adjacency + * flips tile parity, a tile's halo is never written by another thread in the same + * pass, so the update is in-place, race-free, and (being a consistent ordering) + * provably convergent for 0 +#include +#include + +/* Internal-header helpers (poisson_solver_compute_3d_bounds, resolve_omega, etc.) + * are C-linkage inline functions; include under extern "C" since this is a .cu TU. */ +extern "C" { +#include "../linear_solver_internal.h" +} + +/* Tile dimensions swept sequentially by each thread. 8x8 keeps the stale-halo + * fraction small enough that the auto optimal-omega remains stable (mirrors the + * AVX2 Block SOR, whose only staleness is the intra-block left neighbor). */ +#define SOR_GPU_TILE_W 8 +#define SOR_GPU_TILE_H 8 + +/* ============================================================================ + * CONTEXT + * ============================================================================ */ + +typedef struct { + size_t nx, ny, nz, size; + size_t stride_z; + int k_start, k_end; + double inv_dx2, inv_dy2, inv_dz2; + double factor, inv_factor; + double omega; + int block_x, block_y; + int tile_w, tile_h; + + double* d_x; /* solution (in/out) — updated in place */ + double* d_rhs; /* right-hand side */ + double* d_scalar; /* single-double reduction accumulator */ + cudaStream_t stream; +} poisson_sor_gpu_ctx; + +/* ============================================================================ + * HELPERS + * ============================================================================ */ + +/* L2 norm of the Poisson residual rhs - A*x for the field d_field. Returns -1.0 + * on any CUDA failure so the caller can fall back to the fixed iteration cap. */ +static double sor_residual_norm(poisson_sor_gpu_ctx* ctx, const double* d_field, + dim3 cell_grid, dim3 block) { + size_t shmem = (size_t)block.x * block.y * sizeof(double); + if (cudaMemsetAsync(ctx->d_scalar, 0, sizeof(double), ctx->stream) != cudaSuccess) + return -1.0; + lin_gpu_kernel_residual_sq<<stream>>>( + d_field, ctx->d_rhs, ctx->d_scalar, ctx->nx, ctx->ny, + ctx->stride_z, ctx->k_start, ctx->k_end, + ctx->inv_dx2, ctx->inv_dy2, ctx->inv_dz2, ctx->factor); + if (cudaGetLastError() != cudaSuccess) + return -1.0; + double h_sumsq = 0.0; + if (cudaMemcpyAsync(&h_sumsq, ctx->d_scalar, sizeof(double), + cudaMemcpyDeviceToHost, ctx->stream) != cudaSuccess) + return -1.0; + if (cudaStreamSynchronize(ctx->stream) != cudaSuccess) + return -1.0; + return sqrt(h_sumsq); +} + +/* ============================================================================ + * INTERFACE IMPLEMENTATION + * ============================================================================ */ + +static cfd_status_t sor_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) { + if (!gpu_is_available()) { + cfd_set_error(CFD_ERROR_UNSUPPORTED, "CUDA GPU not available at runtime"); + return CFD_ERROR_UNSUPPORTED; + } + + poisson_sor_gpu_ctx* ctx = + (poisson_sor_gpu_ctx*)calloc(1, sizeof(poisson_sor_gpu_ctx)); + if (!ctx) { + cfd_set_error(CFD_ERROR_NOMEM, "Failed to allocate GPU SOR solver context"); + return CFD_ERROR_NOMEM; + } + + ctx->nx = nx; + ctx->ny = ny; + ctx->nz = nz; + ctx->size = nx * ny * (nz > 0 ? nz : 1); + ctx->inv_dx2 = 1.0 / (dx * dx); + ctx->inv_dy2 = 1.0 / (dy * dy); + ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz); + ctx->factor = 2.0 * (ctx->inv_dx2 + ctx->inv_dy2 + ctx->inv_dz2); + ctx->inv_factor = 1.0 / ctx->factor; + ctx->omega = poisson_solver_resolve_omega( + params ? params->omega : 0.0, nx, ny, nz, dx, dy, dz); + + size_t sz, ks, ke; + poisson_solver_compute_3d_bounds(nz, nx, ny, &sz, &ks, &ke); + ctx->stride_z = sz; + ctx->k_start = (int)ks; + ctx->k_end = (int)ke - 1; /* primitives use inclusive k_end */ + + gpu_config_t cfg = gpu_config_default(); + ctx->block_x = cfg.block_size_x; + ctx->block_y = cfg.block_size_y; + ctx->tile_w = SOR_GPU_TILE_W; + ctx->tile_h = SOR_GPU_TILE_H; + + size_t bytes = ctx->size * sizeof(double); + bool ok = cudaMalloc(&ctx->d_x, bytes) == cudaSuccess + && cudaMalloc(&ctx->d_rhs, bytes) == cudaSuccess + && cudaMalloc(&ctx->d_scalar, sizeof(double)) == cudaSuccess + && cudaStreamCreate(&ctx->stream) == cudaSuccess; + if (!ok) { + cudaFree(ctx->d_x); + cudaFree(ctx->d_rhs); + cudaFree(ctx->d_scalar); + free(ctx); + cfd_set_error(CFD_ERROR_NOMEM, "GPU SOR: device allocation failed"); + return CFD_ERROR_NOMEM; + } + + solver->context = ctx; + return CFD_SUCCESS; +} + +static void sor_gpu_destroy(poisson_solver_t* solver) { + if (!solver || !solver->context) + return; + poisson_sor_gpu_ctx* ctx = (poisson_sor_gpu_ctx*)solver->context; + if (ctx->stream) + cudaStreamDestroy(ctx->stream); + cudaFree(ctx->d_x); + cudaFree(ctx->d_rhs); + cudaFree(ctx->d_scalar); + free(ctx); + solver->context = NULL; +} + +static cfd_status_t sor_gpu_solve(poisson_solver_t* solver, + double* x, double* x_temp, const double* rhs, + poisson_solver_stats_t* stats) { + (void)x_temp; /* host scratch unused: the solve runs in-place on device */ + if (!solver || !x || !rhs) + return CFD_ERROR_INVALID; + poisson_sor_gpu_ctx* ctx = (poisson_sor_gpu_ctx*)solver->context; + if (!ctx) + return CFD_ERROR_INVALID; + + const poisson_solver_params_t* p = &solver->params; + size_t nx = ctx->nx, ny = ctx->ny, nz = ctx->nz; + size_t bytes = ctx->size * sizeof(double); + + dim3 block((unsigned)ctx->block_x, (unsigned)ctx->block_y); + /* Sweep grid: one thread per tile. Residual/BC grid: one thread per cell. */ + unsigned n_tiles_x = ((unsigned)(nx - 2) + (unsigned)ctx->tile_w - 1) / (unsigned)ctx->tile_w; + unsigned n_tiles_y = ((unsigned)(ny - 2) + (unsigned)ctx->tile_h - 1) / (unsigned)ctx->tile_h; + dim3 sweep_grid((n_tiles_x + block.x - 1) / block.x, + (n_tiles_y + block.y - 1) / block.y); + dim3 cell_grid((unsigned)((nx - 2 + block.x - 1) / block.x), + (unsigned)((ny - 2 + block.y - 1) / block.y)); + + if (cudaMemcpyAsync(ctx->d_x, x, bytes, cudaMemcpyHostToDevice, ctx->stream) != cudaSuccess + || cudaMemcpyAsync(ctx->d_rhs, rhs, bytes, cudaMemcpyHostToDevice, ctx->stream) + != cudaSuccess) { + return CFD_ERROR; + } + + /* Enforce the Neumann BC on the (possibly warm-started) initial guess before + * measuring r0: the residual stencil reads boundary-adjacent cells. */ + bc_apply_scalar_3d_gpu(ctx->d_x, nx, ny, nz, BC_TYPE_NEUMANN, ctx->stream); + double r0 = sor_residual_norm(ctx, ctx->d_x, cell_grid, block); + + const double RES_FLOOR = 1e-30; + double tol_abs = p->absolute_tolerance; + int can_check = std::isfinite(r0) && (r0 >= 0.0); + double tol_target = can_check ? p->tolerance * r0 : 0.0; + if (tol_target < tol_abs) + tol_target = tol_abs; + + int max_iter = p->max_iterations; + /* Each convergence check is a device-side reduction + stream sync, so never + * poll more often than every CHECK_FLOOR iterations; honor a larger explicit + * check_interval. A small max_iter still polls once at the final iteration. */ + const int CHECK_FLOOR = 20; + int check_every = p->check_interval > CHECK_FLOOR ? p->check_interval : CHECK_FLOOR; + if (max_iter > 0 && check_every > max_iter) + check_every = max_iter; + + double res = r0; + int iter = 0; + int converged = 0; + + if (can_check && r0 <= tol_abs) { + converged = 1; /* already converged */ + } else { + for (iter = 0; iter < max_iter; iter++) { + /* Red tile pass, then black tile pass: the launch boundary is the color + * sync. In-place, so no double buffer. BCs applied after both passes, + * matching the CPU/Red-Black reference. */ + lin_gpu_kernel_block_sor_tile_sweep<<stream>>>( + ctx->d_x, ctx->d_rhs, /*color=*/0, ctx->omega, ctx->tile_w, ctx->tile_h, + nx, ny, ctx->stride_z, ctx->k_start, ctx->k_end, + ctx->inv_dx2, ctx->inv_dy2, ctx->inv_dz2, ctx->inv_factor); + lin_gpu_kernel_block_sor_tile_sweep<<stream>>>( + ctx->d_x, ctx->d_rhs, /*color=*/1, ctx->omega, ctx->tile_w, ctx->tile_h, + nx, ny, ctx->stride_z, ctx->k_start, ctx->k_end, + ctx->inv_dx2, ctx->inv_dy2, ctx->inv_dz2, ctx->inv_factor); + bc_apply_scalar_3d_gpu(ctx->d_x, nx, ny, nz, BC_TYPE_NEUMANN, ctx->stream); + + if (can_check && (iter + 1) % check_every == 0) { + double rnorm = sor_residual_norm(ctx, ctx->d_x, cell_grid, block); + if (rnorm < 0.0 || !std::isfinite(rnorm)) { + can_check = 0; /* residual eval failed: run out the cap */ + } else { + res = rnorm; + if (rnorm <= tol_target || rnorm <= RES_FLOOR) { + converged = 1; + iter++; + break; + } + } + } + } + } + + if (cudaMemcpyAsync(x, ctx->d_x, bytes, cudaMemcpyDeviceToHost, ctx->stream) != cudaSuccess) + return CFD_ERROR; + if (cudaStreamSynchronize(ctx->stream) != cudaSuccess) + return CFD_ERROR; + + if (stats) { + /* sor_residual_norm() returns -1.0 when the device reduction fails; + * normalize such failures to NaN so stats never report a nonsensical + * negative residual. */ + stats->initial_residual = (std::isfinite(r0) && r0 >= 0.0) ? r0 : NAN; + stats->final_residual = (std::isfinite(res) && res >= 0.0) ? res : NAN; + stats->iterations = iter; + stats->status = converged ? POISSON_CONVERGED : POISSON_MAX_ITER; + } + return converged ? CFD_SUCCESS : CFD_ERROR_MAX_ITER; +} + +/* ============================================================================ + * FACTORY + * ============================================================================ */ + +extern "C" poisson_solver_t* create_sor_gpu_solver(void) { + poisson_solver_t* solver = (poisson_solver_t*)calloc(1, sizeof(poisson_solver_t)); + if (!solver) { + cfd_set_error(CFD_ERROR_NOMEM, "Failed to allocate GPU SOR solver"); + return NULL; + } + + solver->name = "sor_gpu"; + solver->description = "SOR iteration (CUDA GPU, Block SOR)"; + solver->method = POISSON_METHOD_SOR; + solver->backend = POISSON_BACKEND_GPU; + solver->params = poisson_solver_params_default(); + + solver->init = sor_gpu_init; + solver->destroy = sor_gpu_destroy; + solver->solve = sor_gpu_solve; + solver->iterate = NULL; /* full solve is implemented directly */ + solver->apply_bc = NULL; /* Neumann applied internally on-device */ + + return solver; +} diff --git a/lib/src/solvers/linear/linear_solver.c b/lib/src/solvers/linear/linear_solver.c index d5744102..6dc67174 100644 --- a/lib/src/solvers/linear/linear_solver.c +++ b/lib/src/solvers/linear/linear_solver.c @@ -177,10 +177,14 @@ poisson_solver_t* poisson_solver_create( switch (backend) { case POISSON_BACKEND_SIMD: return create_sor_simd_solver(); +#ifdef CFD_HAS_CUDA + case POISSON_BACKEND_GPU: + return create_sor_gpu_solver(); +#endif case POISSON_BACKEND_SCALAR: return create_sor_scalar_solver(); default: - return NULL; /* SOR not available for OMP/GPU */ + return NULL; /* SOR not available for OMP */ } case POISSON_METHOD_REDBLACK_SOR: diff --git a/lib/src/solvers/linear/linear_solver_internal.h b/lib/src/solvers/linear/linear_solver_internal.h index eb733668..b78801ab 100644 --- a/lib/src/solvers/linear/linear_solver_internal.h +++ b/lib/src/solvers/linear/linear_solver_internal.h @@ -55,6 +55,7 @@ poisson_solver_t* create_jacobi_gpu_solver(void); poisson_solver_t* create_cg_gpu_solver(void); poisson_solver_t* create_bicgstab_gpu_solver(void); poisson_solver_t* create_redblack_gpu_solver(void); +poisson_solver_t* create_sor_gpu_solver(void); #endif /* BiCGSTAB solvers (for non-symmetric systems) */ diff --git a/tests/math/test_poisson_sor_gpu.c b/tests/math/test_poisson_sor_gpu.c new file mode 100644 index 00000000..5606dab9 --- /dev/null +++ b/tests/math/test_poisson_sor_gpu.c @@ -0,0 +1,268 @@ +/** + * @file test_poisson_sor_gpu.c + * @brief GPU plain SOR Poisson solver: convergence + consistency vs CPU + * + * Verifies the standalone CUDA plain-SOR backend (POISSON_METHOD_SOR + + * POISSON_BACKEND_GPU) that plugs into the poisson_solver_t interface. The GPU + * backend is a "Block SOR": each thread sweeps a small tile sequentially, with + * red-black *tile* coloring (red pass then black pass per iteration), so a tile's + * halo is read from the opposite-color pass rather than written concurrently. It + * is therefore NOT a + * bit-for-bit reproduction of the CPU lexicographic sweep — but it solves the + * same discrete linear system, so once converged it lands on the same field + * (up to the additive Neumann constant) as the CPU reference. + * + * Manufactured solution: p = cos(pi x) cos(pi y) on [0,1]^2. + * nabla^2 p = -2 pi^2 cos(pi x) cos(pi y) -> RHS + * dp/dn = 0 on every face -> exactly satisfies the solver's native Neumann BC. + * The pure-Neumann solution is unique only up to an additive constant, so all + * field comparisons are made after subtracting the interior mean. + * + * The test skips gracefully (no failure) when CUDA is not compiled in (create + * returns NULL) or no GPU device is present at runtime (init returns + * CFD_ERROR_UNSUPPORTED), per the project's optional-backend testing policy. + */ + +#include "unity.h" +#include "cfd/solvers/poisson_solver.h" +#include "cfd/core/memory.h" +#include "cfd/core/indexing.h" + +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define DOMAIN_MIN 0.0 +#define DOMAIN_MAX 1.0 + +void setUp(void) {} +void tearDown(void) {} + +/* ---- helpers ------------------------------------------------------------- */ + +static double* create_field(size_t n) { + return (double*)cfd_calloc(n, sizeof(double)); +} + +static double manufactured_p(double x, double y) { + return cos(M_PI * x) * cos(M_PI * y); +} + +static double manufactured_rhs(double x, double y) { + return -2.0 * M_PI * M_PI * cos(M_PI * x) * cos(M_PI * y); +} + +static void init_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)] = manufactured_rhs(x, y); + } + } +} + +/* Mean of the interior points (used to remove the additive Neumann constant). */ +static double interior_mean(const double* f, size_t nx, size_t ny) { + double sum = 0.0; + size_t count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + sum += f[IDX_2D(i, j, nx)]; + count++; + } + } + return count ? sum / (double)count : 0.0; +} + +/* Max abs interior difference between two fields after removing each one's mean. */ +static double max_diff_demeaned(const double* a, const double* b, size_t nx, size_t ny) { + double ma = interior_mean(a, nx, ny); + double mb = interior_mean(b, nx, ny); + double max_d = 0.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 d = fabs((a[idx] - ma) - (b[idx] - mb)); + if (d > max_d) max_d = d; + } + } + return max_d; +} + +/* L2 error vs the analytical manufactured solution (after mean removal). */ +static double l2_error_vs_analytical(const double* p, size_t nx, size_t ny, + double dx, double dy) { + double mp = interior_mean(p, nx, ny); + double sum = 0.0; + size_t count = 0; + for (size_t j = 1; j < ny - 1; j++) { + for (size_t i = 1; i < nx - 1; i++) { + sum += manufactured_p(DOMAIN_MIN + i * dx, DOMAIN_MIN + j * dy); + count++; + } + } + double ma = count ? sum / (double)count : 0.0; + + double sq = 0.0; + for (size_t j = 1; j < ny - 1; j++) { + double y = DOMAIN_MIN + j * dy; + for (size_t i = 1; i < nx - 1; i++) { + double x = DOMAIN_MIN + i * dx; + size_t idx = IDX_2D(i, j, nx); + double err = (p[idx] - mp) - (manufactured_p(x, y) - ma); + sq += err * err; + } + } + return count ? sqrt(sq / (double)count) : 0.0; +} + +/* Solve the manufactured problem with the given method+backend. Returns: + * 1 on success (solved, *p populated) + * 0 to SKIP (backend unavailable) + * -1 on hard failure (caller should assert). */ +static int solve_backend(poisson_solver_method_t method, poisson_solver_backend_t backend, + double* p, const double* rhs, size_t nx, size_t ny, + double dx, double dy, poisson_solver_stats_t* stats) { + poisson_solver_t* solver = poisson_solver_create(method, backend); + if (!solver) { + return 0; /* backend not compiled in */ + } + + poisson_solver_params_t params = poisson_solver_params_default(); + params.tolerance = 1e-7; + params.absolute_tolerance = 1e-12; + params.max_iterations = 30000; /* generous: Block SOR halos slow convergence */ + + cfd_status_t st = poisson_solver_init(solver, nx, ny, 1, dx, dy, 0.0, ¶ms); + if (st == CFD_ERROR_UNSUPPORTED) { + poisson_solver_destroy(solver); + return 0; /* no GPU device at runtime */ + } + if (st != CFD_SUCCESS) { + poisson_solver_destroy(solver); + return -1; + } + + double* p_temp = create_field(nx * ny); + if (!p_temp) { + poisson_solver_destroy(solver); + return -1; + } + + cfd_status_t solve_st = poisson_solver_solve(solver, p, p_temp, rhs, stats); + cfd_free(p_temp); + poisson_solver_destroy(solver); + + /* CFD_SUCCESS or CFD_ERROR_MAX_ITER (residual still good) both acceptable. */ + return (solve_st == CFD_SUCCESS || solve_st == CFD_ERROR_MAX_ITER) ? 1 : -1; +} + +/* ---- tests --------------------------------------------------------------- */ + +/* GPU Block SOR converges with the auto-resolved optimal omega and drives the + * residual down by many orders of magnitude. This is the key stability check: + * the 2D-tile scheme keeps the stale-halo fraction small enough that omega ~1.8 + * does not diverge (a one-thread-per-row scheme would). */ +void test_sor_gpu_converges(void) { + printf("\n GPU plain SOR: convergence with auto omega...\n"); + size_t nx = 33, ny = 33; + double dx = (DOMAIN_MAX - DOMAIN_MIN) / (nx - 1); + double dy = (DOMAIN_MAX - DOMAIN_MIN) / (ny - 1); + + double* p_gpu = create_field(nx * ny); + double* rhs = create_field(nx * ny); + TEST_ASSERT_NOT_NULL(p_gpu); + TEST_ASSERT_NOT_NULL(rhs); + init_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_stats_t sg = poisson_solver_stats_default(); + int rc = solve_backend(POISSON_METHOD_SOR, POISSON_BACKEND_GPU, + p_gpu, rhs, nx, ny, dx, dy, &sg); + if (rc == 0) { + printf(" SKIPPED (GPU backend unavailable)\n"); + cfd_free(p_gpu); + cfd_free(rhs); + return; + } + TEST_ASSERT_EQUAL_INT_MESSAGE(1, rc, "GPU SOR solve failed"); + printf(" iters=%d init_res=%.3e final_res=%.3e status=%d\n", + sg.iterations, sg.initial_residual, sg.final_residual, sg.status); + + TEST_ASSERT_EQUAL_INT_MESSAGE(POISSON_CONVERGED, sg.status, + "GPU Block SOR did not converge with auto omega (possible instability)"); + TEST_ASSERT_TRUE_MESSAGE(sg.final_residual < 1e-3 * sg.initial_residual, + "GPU SOR did not reduce the residual significantly"); + + cfd_free(p_gpu); + cfd_free(rhs); +} + +/* GPU Block SOR reaches the same discretization-floor accuracy as the CPU plain + * SOR reference. Both solve the identical discrete system, so on convergence they + * agree on the field (up to the additive Neumann constant). The block-staleness + * difference is bounded by the converged tolerance, so a loose ceiling applies. */ +void test_sor_gpu_matches_cpu(void) { + printf("\n GPU plain SOR vs CPU plain SOR consistency...\n"); + size_t nx = 33, ny = 33; + double dx = (DOMAIN_MAX - DOMAIN_MIN) / (nx - 1); + double dy = (DOMAIN_MAX - DOMAIN_MIN) / (ny - 1); + + double* rhs = create_field(nx * ny); + double* p_gpu = create_field(nx * ny); + double* p_cpu = create_field(nx * ny); + TEST_ASSERT_NOT_NULL(rhs); + TEST_ASSERT_NOT_NULL(p_gpu); + TEST_ASSERT_NOT_NULL(p_cpu); + init_rhs(rhs, nx, ny, dx, dy); + + poisson_solver_stats_t sg = poisson_solver_stats_default(); + int rc_gpu = solve_backend(POISSON_METHOD_SOR, POISSON_BACKEND_GPU, + p_gpu, rhs, nx, ny, dx, dy, &sg); + if (rc_gpu == 0) { + printf(" SKIPPED (GPU backend unavailable)\n"); + cfd_free(rhs); + cfd_free(p_gpu); + cfd_free(p_cpu); + return; + } + TEST_ASSERT_EQUAL_INT_MESSAGE(1, rc_gpu, "GPU SOR solve failed"); + + poisson_solver_stats_t sc = poisson_solver_stats_default(); + int rc_cpu = solve_backend(POISSON_METHOD_SOR, POISSON_BACKEND_SCALAR, + p_cpu, rhs, nx, ny, dx, dy, &sc); + TEST_ASSERT_EQUAL_INT_MESSAGE(1, rc_cpu, "CPU SOR reference solve failed"); + + double diff = max_diff_demeaned(p_gpu, p_cpu, nx, ny); + double l2_gpu = l2_error_vs_analytical(p_gpu, nx, ny, dx, dy); + double l2_cpu = l2_error_vs_analytical(p_cpu, nx, ny, dx, dy); + printf(" GPU iters=%d CPU iters=%d max|GPU-CPU|(demeaned)=%.3e\n", + sg.iterations, sc.iterations, diff); + printf(" L2_gpu=%.3e L2_cpu=%.3e\n", l2_gpu, l2_cpu); + + /* Same discrete system: converged fields agree up to the Neumann constant. + * Loose bounds absorb the block-staleness path difference at solver tol. */ + TEST_ASSERT_TRUE_MESSAGE(diff < 1e-3, + "GPU and CPU SOR solutions diverge beyond 1e-3"); + TEST_ASSERT_TRUE_MESSAGE(fabs(l2_gpu - l2_cpu) < 1e-3, + "GPU SOR accuracy differs from CPU reference"); + TEST_ASSERT_TRUE_MESSAGE(l2_gpu < 1e-1, + "GPU SOR L2 error vs analytical unexpectedly large"); + + cfd_free(rhs); + cfd_free(p_gpu); + cfd_free(p_cpu); +} + +int main(void) { + UNITY_BEGIN(); + printf("\n========================================\n"); + printf("GPU PLAIN SOR (POISSON) SOLVER TESTS\n"); + printf("========================================\n"); + RUN_TEST(test_sor_gpu_converges); + RUN_TEST(test_sor_gpu_matches_cpu); + return UNITY_END(); +}