Line data    Source code

       1              : /**
       2              :  * @file linear_solver_jacobi.c
       3              :  * @brief Jacobi iteration solver - scalar CPU implementation
       4              :  *
       5              :  * Jacobi method characteristics:
       6              :  * - Fully parallelizable (reads from old array, writes to new array)
       7              :  * - Requires temporary buffer for double-buffering
       8              :  * - Slower convergence than SOR (~2x iterations)
       9              :  * - Better suited for GPU/massively parallel execution
      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              : #include <string.h>
      20              : 
      21              : /* ============================================================================
      22              :  * JACOBI CONTEXT
      23              :  * ============================================================================ */
      24              : 
      25              : typedef struct {
      26              :     double dx2;        /* dx^2 */
      27              :     double dy2;        /* dy^2 */
      28              :     double inv_dz2;
      29              :     double inv_factor; /* 1 / (2 * (1/dx^2 + 1/dy^2)) */
      30              :     size_t stride_z;
      31              :     size_t k_start;
      32              :     size_t k_end;
      33              :     size_t nz;
      34              :     int initialized;
      35              : } jacobi_context_t;
      36              : 
      37              : /* ============================================================================
      38              :  * JACOBI SCALAR IMPLEMENTATION
      39              :  * ============================================================================ */
      40              : 
      41           25 : static cfd_status_t jacobi_scalar_init(
      42              :     poisson_solver_t* solver,
      43              :     size_t nx, size_t ny, size_t nz,
      44              :     double dx, double dy, double dz,
      45              :     const poisson_solver_params_t* params)
      46              : {
      47           25 :     (void)params;
      48              : 
      49           25 :     jacobi_context_t* ctx = (jacobi_context_t*)cfd_calloc(1, sizeof(jacobi_context_t));
      50           25 :     if (!ctx) {
      51              :         return CFD_ERROR_NOMEM;
      52              :     }
      53              : 
      54           25 :     ctx->dx2 = dx * dx;
      55           25 :     ctx->dy2 = dy * dy;
      56           25 :     ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
      57           25 :     ctx->nz = nz;
      58           25 :     poisson_solver_compute_3d_bounds(nz, nx, ny,
      59              :         &ctx->stride_z, &ctx->k_start, &ctx->k_end);
      60              : 
      61           25 :     double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
      62           25 :     ctx->inv_factor = 1.0 / factor;
      63           25 :     ctx->initialized = 1;
      64              : 
      65           25 :     solver->context = ctx;
      66           25 :     return CFD_SUCCESS;
      67              : }
      68              : 
      69           26 : static void jacobi_scalar_destroy(poisson_solver_t* solver) {
      70           26 :     if (solver && solver->context) {
      71           25 :         cfd_free(solver->context);
      72           25 :         solver->context = NULL;
      73              :     }
      74           26 : }
      75              : 
      76        70696 : static cfd_status_t jacobi_scalar_iterate(
      77              :     poisson_solver_t* solver,
      78              :     double* x,
      79              :     double* x_temp,
      80              :     const double* rhs,
      81              :     double* residual)
      82              : {
      83        70696 :     if (!x_temp) {
      84              :         return CFD_ERROR_INVALID;  /* Jacobi requires temp buffer */
      85              :     }
      86              : 
      87        70696 :     jacobi_context_t* ctx = (jacobi_context_t*)solver->context;
      88        70696 :     size_t nx = solver->nx;
      89        70696 :     size_t ny = solver->ny;
      90        70696 :     double dx2 = ctx->dx2;
      91        70696 :     double dy2 = ctx->dy2;
      92        70696 :     double inv_factor = ctx->inv_factor;
      93              : 
      94        70696 :     double* p_old = x;
      95        70696 :     double* p_new = x_temp;
      96              : 
      97        70696 :     size_t stride_z = ctx->stride_z;
      98        70696 :     double inv_dz2 = ctx->inv_dz2;
      99              : 
     100              :     /* Jacobi update: reads from p_old, writes to p_new */
     101       167992 :     for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
     102      2814072 :         for (size_t j = 1; j < ny - 1; j++) {
     103     97802832 :             for (size_t i = 1; i < nx - 1; i++) {
     104     95086056 :                 size_t idx = k * stride_z + IDX_2D(i, j, nx);
     105              : 
     106     95086056 :                 double p_result = -(rhs[idx]
     107     95086056 :                     - (p_old[idx + 1] + p_old[idx - 1]) / dx2
     108     95086056 :                     - (p_old[idx + nx] + p_old[idx - nx]) / dy2
     109     95086056 :                     - (p_old[idx + stride_z] + p_old[idx - stride_z]) * inv_dz2
     110              :                     ) * inv_factor;
     111              : 
     112     95086056 :                 p_new[idx] = p_result;
     113              :             }
     114              :         }
     115              :     }
     116              : 
     117              :     /* Copy result back to x */
     118        70696 :     memcpy(x, x_temp, nx * ny * ctx->nz * sizeof(double));
     119              : 
     120              :     /* Apply boundary conditions */
     121        70696 :     poisson_solver_apply_bc(solver, x);
     122              : 
     123              :     /* Compute residual if requested */
     124        70696 :     if (residual) {
     125        70696 :         *residual = poisson_solver_compute_residual(solver, x, rhs);
     126              :     }
     127              : 
     128              :     return CFD_SUCCESS;
     129              : }
     130              : 
     131              : /* ============================================================================
     132              :  * FACTORY FUNCTION
     133              :  * ============================================================================ */
     134              : 
     135           26 : poisson_solver_t* create_jacobi_scalar_solver(void) {
     136           26 :     poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
     137           26 :     if (!solver) {
     138              :         return NULL;
     139              :     }
     140              : 
     141           26 :     solver->name = POISSON_SOLVER_TYPE_JACOBI_SCALAR;
     142           26 :     solver->description = "Jacobi iteration (scalar CPU)";
     143           26 :     solver->method = POISSON_METHOD_JACOBI;
     144           26 :     solver->backend = POISSON_BACKEND_SCALAR;
     145           26 :     solver->params = poisson_solver_params_default();
     146           26 :     solver->params.max_iterations = 2000;  /* Jacobi needs more iterations */
     147           26 :     solver->params.check_interval = 10;    /* Check less frequently for speed */
     148              : 
     149           26 :     solver->init = jacobi_scalar_init;
     150           26 :     solver->destroy = jacobi_scalar_destroy;
     151           26 :     solver->solve = NULL;  /* Use common solve loop */
     152           26 :     solver->iterate = jacobi_scalar_iterate;
     153           26 :     solver->apply_bc = NULL;  /* Use default Neumann */
     154              : 
     155           26 :     return solver;
     156              : }