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           21 : 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           21 :     sor_context_t* ctx = (sor_context_t*)cfd_calloc(1, sizeof(sor_context_t));
      47           21 :     if (!ctx) {
      48              :         return CFD_ERROR_NOMEM;
      49              :     }
      50              : 
      51           21 :     ctx->dx2 = dx * dx;
      52           21 :     ctx->dy2 = dy * dy;
      53           21 :     ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
      54           21 :     poisson_solver_compute_3d_bounds(nz, nx, ny,
      55              :         &ctx->stride_z, &ctx->k_start, &ctx->k_end);
      56              : 
      57           21 :     double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
      58           21 :     ctx->inv_factor = 1.0 / factor;
      59           21 :     ctx->omega = params ? params->omega : 1.5;
      60           21 :     ctx->initialized = 1;
      61              : 
      62           21 :     solver->context = ctx;
      63           21 :     return CFD_SUCCESS;
      64              : }
      65              : 
      66           22 : static void sor_scalar_destroy(poisson_solver_t* solver) {
      67           22 :     if (solver && solver->context) {
      68           21 :         cfd_free(solver->context);
      69           21 :         solver->context = NULL;
      70              :     }
      71           22 : }
      72              : 
      73              : /**
      74              :  * Sequential SOR iteration
      75              :  *
      76              :  * Updates cells in row-major order, using already-updated values
      77              :  * from the current iteration (Gauss-Seidel behavior).
      78              :  */
      79        57021 : static cfd_status_t sor_scalar_iterate(
      80              :     poisson_solver_t* solver,
      81              :     double* x,
      82              :     double* x_temp,
      83              :     const double* rhs,
      84              :     double* residual)
      85              : {
      86        57021 :     (void)x_temp;  /* Not needed for in-place SOR */
      87              : 
      88        57021 :     sor_context_t* ctx = (sor_context_t*)solver->context;
      89        57021 :     size_t nx = solver->nx;
      90        57021 :     size_t ny = solver->ny;
      91        57021 :     double dx2 = ctx->dx2;
      92        57021 :     double dy2 = ctx->dy2;
      93        57021 :     double inv_factor = ctx->inv_factor;
      94        57021 :     double omega = ctx->omega;
      95              : 
      96        57021 :     size_t stride_z = ctx->stride_z;
      97        57021 :     double inv_dz2 = ctx->inv_dz2;
      98              : 
      99              :     /* Single sweep: row-major order */
     100       118214 :     for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
     101      2046456 :         for (size_t j = 1; j < ny - 1; j++) {
     102     80209616 :             for (size_t i = 1; i < nx - 1; i++) {
     103     78224353 :                 size_t idx = k * stride_z + IDX_2D(i, j, nx);
     104              : 
     105              :                 /* Compute Gauss-Seidel update
     106              :                  * Note: x[idx-1] and x[idx-nx] are already updated this iteration
     107              :                  */
     108     78224353 :                 double p_new = -(rhs[idx]
     109     78224353 :                     - (x[idx + 1] + x[idx - 1]) / dx2
     110     78224353 :                     - (x[idx + nx] + x[idx - nx]) / dy2
     111     78224353 :                     - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
     112              :                     ) * inv_factor;
     113              : 
     114              :                 /* SOR relaxation */
     115     78224353 :                 x[idx] = x[idx] + omega * (p_new - x[idx]);
     116              :             }
     117              :         }
     118              :     }
     119              : 
     120              :     /* Apply boundary conditions */
     121        57021 :     poisson_solver_apply_bc(solver, x);
     122              : 
     123              :     /* Compute residual if requested */
     124        57021 :     if (residual) {
     125        57021 :         *residual = poisson_solver_compute_residual(solver, x, rhs);
     126              :     }
     127              : 
     128        57021 :     return CFD_SUCCESS;
     129              : }
     130              : 
     131              : /* ============================================================================
     132              :  * FACTORY FUNCTION
     133              :  * ============================================================================ */
     134              : 
     135           22 : poisson_solver_t* create_sor_scalar_solver(void) {
     136           22 :     poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
     137           22 :     if (!solver) {
     138              :         return NULL;
     139              :     }
     140              : 
     141           22 :     solver->name = POISSON_SOLVER_TYPE_SOR_SCALAR;
     142           22 :     solver->description = "SOR iteration (scalar CPU, sequential)";
     143           22 :     solver->method = POISSON_METHOD_SOR;
     144           22 :     solver->backend = POISSON_BACKEND_SCALAR;
     145           22 :     solver->params = poisson_solver_params_default();
     146              : 
     147           22 :     solver->init = sor_scalar_init;
     148           22 :     solver->destroy = sor_scalar_destroy;
     149           22 :     solver->solve = NULL;  /* Use common solve loop */
     150           22 :     solver->iterate = sor_scalar_iterate;
     151           22 :     solver->apply_bc = NULL;  /* Use default Neumann */
     152              : 
     153           22 :     return solver;
     154              : }