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 12 : 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 12 : redblack_context_t* ctx = (redblack_context_t*)cfd_calloc(1, sizeof(redblack_context_t));
47 12 : if (!ctx) {
48 : return CFD_ERROR_NOMEM;
49 : }
50 :
51 12 : ctx->dx2 = dx * dx;
52 12 : ctx->dy2 = dy * dy;
53 12 : ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
54 12 : poisson_solver_compute_3d_bounds(nz, nx, ny,
55 : &ctx->stride_z, &ctx->k_start, &ctx->k_end);
56 :
57 12 : double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
58 12 : ctx->inv_factor = 1.0 / factor;
59 12 : ctx->omega = params ? params->omega : 1.5;
60 12 : ctx->initialized = 1;
61 :
62 12 : solver->context = ctx;
63 12 : return CFD_SUCCESS;
64 : }
65 :
66 14 : static void redblack_scalar_destroy(poisson_solver_t* solver) {
67 14 : if (solver && solver->context) {
68 12 : cfd_free(solver->context);
69 12 : solver->context = NULL;
70 : }
71 14 : }
72 :
73 : /**
74 : * Red-Black SOR iteration (scalar)
75 : *
76 : * Red cells: (i+j) % 2 == 0
77 : * Black cells: (i+j) % 2 == 1
78 : */
79 55662 : static cfd_status_t redblack_scalar_iterate(
80 : poisson_solver_t* solver,
81 : double* x,
82 : double* x_temp,
83 : const double* rhs,
84 : double* residual)
85 : {
86 55662 : (void)x_temp; /* Not needed for in-place SOR */
87 :
88 55662 : redblack_context_t* ctx = (redblack_context_t*)solver->context;
89 55662 : size_t nx = solver->nx;
90 55662 : size_t ny = solver->ny;
91 55662 : double dx2 = ctx->dx2;
92 55662 : double dy2 = ctx->dy2;
93 55662 : double inv_factor = ctx->inv_factor;
94 55662 : double omega = ctx->omega;
95 :
96 55662 : size_t stride_z = ctx->stride_z;
97 55662 : double inv_dz2 = ctx->inv_dz2;
98 :
99 : /* Red sweep: (i+j+k) % 2 == 0 */
100 115720 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
101 2006496 : for (size_t j = 1; j < ny - 1; j++) {
102 1946438 : size_t i_start = ((j + k) % 2 == 0) ? 1 : 2;
103 40403510 : for (size_t i = i_start; i < nx - 1; i += 2) {
104 38457072 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
105 :
106 38457072 : double p_new = -(rhs[idx]
107 38457072 : - (x[idx + 1] + x[idx - 1]) / dx2
108 38457072 : - (x[idx + nx] + x[idx - nx]) / dy2
109 38457072 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
110 : ) * inv_factor;
111 :
112 : /* SOR update */
113 38457072 : x[idx] = x[idx] + omega * (p_new - x[idx]);
114 : }
115 : }
116 : }
117 :
118 : /* Black sweep: (i+j+k) % 2 == 1 */
119 115720 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
120 2006496 : for (size_t j = 1; j < ny - 1; j++) {
121 1946438 : size_t i_start = ((j + k) % 2 == 0) ? 2 : 1;
122 40458544 : for (size_t i = i_start; i < nx - 1; i += 2) {
123 38512106 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
124 :
125 38512106 : double p_new = -(rhs[idx]
126 38512106 : - (x[idx + 1] + x[idx - 1]) / dx2
127 38512106 : - (x[idx + nx] + x[idx - nx]) / dy2
128 38512106 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
129 : ) * inv_factor;
130 :
131 : /* SOR update */
132 38512106 : x[idx] = x[idx] + omega * (p_new - x[idx]);
133 : }
134 : }
135 : }
136 :
137 : /* Apply boundary conditions */
138 55662 : poisson_solver_apply_bc(solver, x);
139 :
140 : /* Compute residual if requested */
141 55662 : if (residual) {
142 55662 : *residual = poisson_solver_compute_residual(solver, x, rhs);
143 : }
144 :
145 55662 : return CFD_SUCCESS;
146 : }
147 :
148 : /* ============================================================================
149 : * FACTORY FUNCTION
150 : * ============================================================================ */
151 :
152 14 : poisson_solver_t* create_redblack_scalar_solver(void) {
153 14 : poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
154 14 : if (!solver) {
155 : return NULL;
156 : }
157 :
158 14 : solver->name = POISSON_SOLVER_TYPE_REDBLACK_SCALAR;
159 14 : solver->description = "Red-Black SOR iteration (scalar CPU)";
160 14 : solver->method = POISSON_METHOD_REDBLACK_SOR;
161 14 : solver->backend = POISSON_BACKEND_SCALAR;
162 14 : solver->params = poisson_solver_params_default();
163 :
164 14 : solver->init = redblack_scalar_init;
165 14 : solver->destroy = redblack_scalar_destroy;
166 14 : solver->solve = NULL; /* Use common solve loop */
167 14 : solver->iterate = redblack_scalar_iterate;
168 14 : solver->apply_bc = NULL;
169 :
170 14 : return solver;
171 : }
|
