Line data    Source code

       1              : /**
       2              :  * @file linear_solver_sor.c
       3              :  * @brief Sequential SOR (Successive Over-Relaxation) solver - scalar CPU implementation
       4              :  *
       5              :  * SOR characteristics:
       6              :  * - Inherently sequential (Gauss-Seidel + relaxation)
       7              :  * - In-place updates (no temporary buffer needed)
       8              :  * - Fastest convergence per iteration
       9              :  * - Cannot be parallelized (unlike Red-Black SOR)
      10              :  */
      11              : 
      12              : #include "../linear_solver_internal.h"
      13              : 
      14              : #include "cfd/boundary/boundary_conditions.h"
      15              : #include "cfd/core/indexing.h"
      16              : #include "cfd/core/memory.h"
      17              : 
      18              : #include <math.h>
      19              : 
      20              : /* ============================================================================
      21              :  * SOR CONTEXT
      22              :  * ============================================================================ */
      23              : 
      24              : typedef struct {
      25              :     double dx2;        /* dx^2 */
      26              :     double dy2;        /* dy^2 */
      27              :     double inv_dz2;
      28              :     double inv_factor; /* 1 / (2 * (1/dx^2 + 1/dy^2)) */
      29              :     double omega;      /* SOR relaxation parameter */
      30              :     size_t stride_z;
      31              :     size_t k_start;
      32              :     size_t k_end;
      33              :     int initialized;
      34              : } sor_context_t;
      35              : 
      36              : /* ============================================================================
      37              :  * SOR SCALAR IMPLEMENTATION
      38              :  * ============================================================================ */
      39              : 
      40           23 : static cfd_status_t sor_scalar_init(
      41              :     poisson_solver_t* solver,
      42              :     size_t nx, size_t ny, size_t nz,
      43              :     double dx, double dy, double dz,
      44              :     const poisson_solver_params_t* params)
      45              : {
      46           23 :     sor_context_t* ctx = (sor_context_t*)cfd_calloc(1, sizeof(sor_context_t));
      47           23 :     if (!ctx) {
      48              :         return CFD_ERROR_NOMEM;
      49              :     }
      50              : 
      51           23 :     ctx->dx2 = dx * dx;
      52           23 :     ctx->dy2 = dy * dy;
      53           23 :     ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
      54           23 :     poisson_solver_compute_3d_bounds(nz, nx, ny,
      55              :         &ctx->stride_z, &ctx->k_start, &ctx->k_end);
      56              : 
      57           23 :     double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
      58           23 :     ctx->inv_factor = 1.0 / factor;
      59           23 :     ctx->omega = poisson_solver_resolve_omega(
      60              :         params ? params->omega : 0.0, nx, ny, nz, dx, dy, dz);
      61           23 :     ctx->initialized = 1;
      62              : 
      63           23 :     solver->context = ctx;
      64           23 :     return CFD_SUCCESS;
      65              : }
      66              : 
      67           24 : static void sor_scalar_destroy(poisson_solver_t* solver) {
      68           24 :     if (solver && solver->context) {
      69           23 :         cfd_free(solver->context);
      70           23 :         solver->context = NULL;
      71              :     }
      72           24 : }
      73              : 
      74              : /**
      75              :  * Sequential SOR iteration
      76              :  *
      77              :  * Updates cells in row-major order, using already-updated values
      78              :  * from the current iteration (Gauss-Seidel behavior).
      79              :  */
      80        57681 : static cfd_status_t sor_scalar_iterate(
      81              :     poisson_solver_t* solver,
      82              :     double* x,
      83              :     double* x_temp,
      84              :     const double* rhs,
      85              :     double* residual)
      86              : {
      87        57681 :     (void)x_temp;  /* Not needed for in-place SOR */
      88              : 
      89        57681 :     sor_context_t* ctx = (sor_context_t*)solver->context;
      90        57681 :     size_t nx = solver->nx;
      91        57681 :     size_t ny = solver->ny;
      92        57681 :     double dx2 = ctx->dx2;
      93        57681 :     double dy2 = ctx->dy2;
      94        57681 :     double inv_factor = ctx->inv_factor;
      95        57681 :     double omega = ctx->omega;
      96              : 
      97        57681 :     size_t stride_z = ctx->stride_z;
      98        57681 :     double inv_dz2 = ctx->inv_dz2;
      99              : 
     100              :     /* Single sweep: row-major order */
     101       119534 :     for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
     102      2057016 :         for (size_t j = 1; j < ny - 1; j++) {
     103     80368016 :             for (size_t i = 1; i < nx - 1; i++) {
     104     78372853 :                 size_t idx = k * stride_z + IDX_2D(i, j, nx);
     105              : 
     106              :                 /* Compute Gauss-Seidel update
     107              :                  * Note: x[idx-1] and x[idx-nx] are already updated this iteration
     108              :                  */
     109     78372853 :                 double p_new = -(rhs[idx]
     110     78372853 :                     - (x[idx + 1] + x[idx - 1]) / dx2
     111     78372853 :                     - (x[idx + nx] + x[idx - nx]) / dy2
     112     78372853 :                     - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
     113              :                     ) * inv_factor;
     114              : 
     115              :                 /* SOR relaxation */
     116     78372853 :                 x[idx] = x[idx] + omega * (p_new - x[idx]);
     117              :             }
     118              :         }
     119              :     }
     120              : 
     121              :     /* Apply boundary conditions */
     122        57681 :     poisson_solver_apply_bc(solver, x);
     123              : 
     124              :     /* Compute residual if requested */
     125        57681 :     if (residual) {
     126        57681 :         *residual = poisson_solver_compute_residual(solver, x, rhs);
     127              :     }
     128              : 
     129        57681 :     return CFD_SUCCESS;
     130              : }
     131              : 
     132              : /* ============================================================================
     133              :  * FACTORY FUNCTION
     134              :  * ============================================================================ */
     135              : 
     136           24 : poisson_solver_t* create_sor_scalar_solver(void) {
     137           24 :     poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
     138           24 :     if (!solver) {
     139              :         return NULL;
     140              :     }
     141              : 
     142           24 :     solver->name = POISSON_SOLVER_TYPE_SOR_SCALAR;
     143           24 :     solver->description = "SOR iteration (scalar CPU, sequential)";
     144           24 :     solver->method = POISSON_METHOD_SOR;
     145           24 :     solver->backend = POISSON_BACKEND_SCALAR;
     146           24 :     solver->params = poisson_solver_params_default();
     147              : 
     148           24 :     solver->init = sor_scalar_init;
     149           24 :     solver->destroy = sor_scalar_destroy;
     150           24 :     solver->solve = NULL;  /* Use common solve loop */
     151           24 :     solver->iterate = sor_scalar_iterate;
     152           24 :     solver->apply_bc = NULL;  /* Use default Neumann */
     153              : 
     154           24 :     return solver;
     155              : }