Line data Source code
1 : /**
2 : * @file linear_solver_redblack.c
3 : * @brief Red-Black SOR solver - scalar CPU implementation
4 : *
5 : * Red-Black SOR characteristics:
6 : * - Two-color sweep (red then black) allows parallelization
7 : * - In-place updates (no temporary buffer needed)
8 : * - SOR acceleration (omega > 1) for faster convergence
9 : * - Best balance of convergence speed and parallelism
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 : * RED-BLACK 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 : } redblack_context_t;
35 :
36 : /* ============================================================================
37 : * RED-BLACK SCALAR IMPLEMENTATION
38 : * ============================================================================ */
39 :
40 20 : static cfd_status_t redblack_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 20 : redblack_context_t* ctx = (redblack_context_t*)cfd_calloc(1, sizeof(redblack_context_t));
47 20 : if (!ctx) {
48 : return CFD_ERROR_NOMEM;
49 : }
50 :
51 20 : ctx->dx2 = dx * dx;
52 20 : ctx->dy2 = dy * dy;
53 20 : ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
54 20 : poisson_solver_compute_3d_bounds(nz, nx, ny,
55 : &ctx->stride_z, &ctx->k_start, &ctx->k_end);
56 :
57 20 : double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
58 20 : ctx->inv_factor = 1.0 / factor;
59 20 : ctx->omega = poisson_solver_resolve_omega(
60 : params ? params->omega : 0.0, nx, ny, nz, dx, dy, dz);
61 20 : ctx->initialized = 1;
62 :
63 20 : solver->context = ctx;
64 20 : return CFD_SUCCESS;
65 : }
66 :
67 22 : static void redblack_scalar_destroy(poisson_solver_t* solver) {
68 22 : if (solver && solver->context) {
69 20 : cfd_free(solver->context);
70 20 : solver->context = NULL;
71 : }
72 22 : }
73 :
74 : /**
75 : * Red-Black SOR iteration (scalar)
76 : *
77 : * Red cells: (i+j) % 2 == 0
78 : * Black cells: (i+j) % 2 == 1
79 : */
80 56732 : static cfd_status_t redblack_scalar_iterate(
81 : poisson_solver_t* solver,
82 : double* x,
83 : double* x_temp,
84 : const double* rhs,
85 : double* residual)
86 : {
87 56732 : (void)x_temp; /* Not needed for in-place SOR */
88 :
89 56732 : redblack_context_t* ctx = (redblack_context_t*)solver->context;
90 56732 : size_t nx = solver->nx;
91 56732 : size_t ny = solver->ny;
92 56732 : double dx2 = ctx->dx2;
93 56732 : double dy2 = ctx->dy2;
94 56732 : double inv_factor = ctx->inv_factor;
95 56732 : double omega = ctx->omega;
96 :
97 56732 : size_t stride_z = ctx->stride_z;
98 56732 : double inv_dz2 = ctx->inv_dz2;
99 :
100 : /* Red sweep: (i+j+k) % 2 == 0 */
101 117860 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
102 2050512 : for (size_t j = 1; j < ny - 1; j++) {
103 1989384 : size_t i_start = ((j + k) % 2 == 0) ? 1 : 2;
104 41293624 : for (size_t i = i_start; i < nx - 1; i += 2) {
105 39304240 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
106 :
107 39304240 : double p_new = -(rhs[idx]
108 39304240 : - (x[idx + 1] + x[idx - 1]) / dx2
109 39304240 : - (x[idx + nx] + x[idx - nx]) / dy2
110 39304240 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
111 : ) * inv_factor;
112 :
113 : /* SOR update */
114 39304240 : x[idx] = x[idx] + omega * (p_new - x[idx]);
115 : }
116 : }
117 : }
118 :
119 : /* Black sweep: (i+j+k) % 2 == 1 */
120 117860 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
121 2050512 : for (size_t j = 1; j < ny - 1; j++) {
122 1989384 : size_t i_start = ((j + k) % 2 == 0) ? 2 : 1;
123 41349728 : for (size_t i = i_start; i < nx - 1; i += 2) {
124 39360344 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
125 :
126 39360344 : double p_new = -(rhs[idx]
127 39360344 : - (x[idx + 1] + x[idx - 1]) / dx2
128 39360344 : - (x[idx + nx] + x[idx - nx]) / dy2
129 39360344 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
130 : ) * inv_factor;
131 :
132 : /* SOR update */
133 39360344 : x[idx] = x[idx] + omega * (p_new - x[idx]);
134 : }
135 : }
136 : }
137 :
138 : /* Apply boundary conditions */
139 56732 : poisson_solver_apply_bc(solver, x);
140 :
141 : /* Compute residual if requested */
142 56732 : if (residual) {
143 56732 : *residual = poisson_solver_compute_residual(solver, x, rhs);
144 : }
145 :
146 56732 : return CFD_SUCCESS;
147 : }
148 :
149 : /* ============================================================================
150 : * FACTORY FUNCTION
151 : * ============================================================================ */
152 :
153 22 : poisson_solver_t* create_redblack_scalar_solver(void) {
154 22 : poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
155 22 : if (!solver) {
156 : return NULL;
157 : }
158 :
159 22 : solver->name = POISSON_SOLVER_TYPE_REDBLACK_SCALAR;
160 22 : solver->description = "Red-Black SOR iteration (scalar CPU)";
161 22 : solver->method = POISSON_METHOD_REDBLACK_SOR;
162 22 : solver->backend = POISSON_BACKEND_SCALAR;
163 22 : solver->params = poisson_solver_params_default();
164 :
165 22 : solver->init = redblack_scalar_init;
166 22 : solver->destroy = redblack_scalar_destroy;
167 22 : solver->solve = NULL; /* Use common solve loop */
168 22 : solver->iterate = redblack_scalar_iterate;
169 22 : solver->apply_bc = NULL;
170 :
171 22 : return solver;
172 : }
|
