Line data Source code
1 : /**
2 : * @file linear_solver_redblack_omp.c
3 : * @brief Red-Black SOR solver - OpenMP parallel 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/core/indexing.h"
15 : #include "cfd/core/memory.h"
16 :
17 : #include <math.h>
18 :
19 : #ifdef CFD_ENABLE_OPENMP
20 : #include <omp.h>
21 : #endif
22 :
23 : /* ============================================================================
24 : * RED-BLACK CONTEXT
25 : * ============================================================================ */
26 :
27 : typedef struct {
28 : double dx2; /* dx^2 */
29 : double dy2; /* dy^2 */
30 : double inv_dz2; /* 1/dz^2 (0 for 2D) */
31 : double inv_factor; /* 1 / (2 * (1/dx^2 + 1/dy^2 + inv_dz2)) */
32 : double omega; /* SOR relaxation parameter */
33 : size_t stride_z; /* nx*ny for 3D, 0 for 2D */
34 : size_t k_start; /* first interior k index */
35 : size_t k_end; /* one-past-last interior k index */
36 : int initialized;
37 : } redblack_omp_context_t;
38 :
39 : /* ============================================================================
40 : * RED-BLACK OPENMP IMPLEMENTATION
41 : * ============================================================================ */
42 :
43 : #ifdef CFD_ENABLE_OPENMP
44 :
45 1 : static cfd_status_t redblack_omp_init(
46 : poisson_solver_t* solver,
47 : size_t nx, size_t ny, size_t nz,
48 : double dx, double dy, double dz,
49 : const poisson_solver_params_t* params)
50 : {
51 1 : redblack_omp_context_t* ctx = (redblack_omp_context_t*)cfd_calloc(1, sizeof(redblack_omp_context_t));
52 1 : if (!ctx) {
53 : return CFD_ERROR_NOMEM;
54 : }
55 :
56 1 : ctx->dx2 = dx * dx;
57 1 : ctx->dy2 = dy * dy;
58 1 : ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
59 1 : poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end);
60 1 : double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
61 1 : ctx->inv_factor = 1.0 / factor;
62 1 : ctx->omega = params ? params->omega : 1.5;
63 1 : ctx->initialized = 1;
64 :
65 1 : solver->context = ctx;
66 1 : return CFD_SUCCESS;
67 : }
68 :
69 2 : static void redblack_omp_destroy(poisson_solver_t* solver) {
70 2 : if (solver && solver->context) {
71 1 : cfd_free(solver->context);
72 1 : solver->context = NULL;
73 : }
74 2 : }
75 :
76 157 : static cfd_status_t redblack_omp_iterate(
77 : poisson_solver_t* solver,
78 : double* x,
79 : double* x_temp,
80 : const double* rhs,
81 : double* residual)
82 : {
83 157 : (void)x_temp;
84 :
85 157 : redblack_omp_context_t* ctx = (redblack_omp_context_t*)solver->context;
86 157 : size_t nx = solver->nx;
87 157 : size_t ny = solver->ny;
88 157 : double dx2 = ctx->dx2;
89 157 : double dy2 = ctx->dy2;
90 157 : double inv_dz2 = ctx->inv_dz2;
91 157 : double inv_factor = ctx->inv_factor;
92 157 : double omega = ctx->omega;
93 157 : size_t stride_z = ctx->stride_z;
94 157 : size_t k_start = ctx->k_start;
95 157 : size_t k_end = ctx->k_end;
96 :
97 : /* Red sweep (parallel) */
98 2512 : for (size_t k = k_start; k < k_end; k++) {
99 2355 : int j;
100 2355 : #pragma omp parallel for schedule(static)
101 : for (j = 1; j < (int)ny - 1; j++) {
102 : size_t i_start = ((j + k) % 2 == 0) ? 1 : 2;
103 : size_t i;
104 : for (i = i_start; i < nx - 1; i += 2) {
105 : size_t idx = k * stride_z + IDX_2D(i, (size_t)j, nx);
106 :
107 : double p_new = -(rhs[idx]
108 : - (x[idx + 1] + x[idx - 1]) / dx2
109 : - (x[idx + nx] + x[idx - nx]) / dy2
110 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
111 : ) * inv_factor;
112 :
113 : x[idx] = x[idx] + omega * (p_new - x[idx]);
114 : }
115 : }
116 : }
117 :
118 : /* Black sweep (parallel) */
119 2512 : for (size_t k = k_start; k < k_end; k++) {
120 2355 : int j;
121 2355 : #pragma omp parallel for schedule(static)
122 : for (j = 1; j < (int)ny - 1; j++) {
123 : size_t i_start = ((j + k) % 2 == 0) ? 2 : 1;
124 : size_t i;
125 : for (i = i_start; i < nx - 1; i += 2) {
126 : size_t idx = k * stride_z + IDX_2D(i, (size_t)j, nx);
127 :
128 : double p_new = -(rhs[idx]
129 : - (x[idx + 1] + x[idx - 1]) / dx2
130 : - (x[idx + nx] + x[idx - nx]) / dy2
131 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
132 : ) * inv_factor;
133 :
134 : x[idx] = x[idx] + omega * (p_new - x[idx]);
135 : }
136 : }
137 : }
138 :
139 : /* Apply boundary conditions */
140 157 : poisson_solver_apply_bc(solver, x);
141 :
142 : /* Compute residual if requested */
143 157 : if (residual) {
144 157 : *residual = poisson_solver_compute_residual(solver, x, rhs);
145 : }
146 :
147 157 : return CFD_SUCCESS;
148 : }
149 :
150 : #endif /* CFD_ENABLE_OPENMP */
151 :
152 : /* ============================================================================
153 : * FACTORY FUNCTION
154 : * ============================================================================ */
155 :
156 : #ifdef CFD_ENABLE_OPENMP
157 2 : poisson_solver_t* create_redblack_omp_solver(void) {
158 2 : poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
159 2 : if (!solver) {
160 : return NULL;
161 : }
162 :
163 2 : solver->name = POISSON_SOLVER_TYPE_REDBLACK_OMP;
164 2 : solver->description = "Red-Black SOR iteration (OpenMP parallel)";
165 2 : solver->method = POISSON_METHOD_REDBLACK_SOR;
166 2 : solver->backend = POISSON_BACKEND_OMP;
167 2 : solver->params = poisson_solver_params_default();
168 :
169 2 : solver->init = redblack_omp_init;
170 2 : solver->destroy = redblack_omp_destroy;
171 2 : solver->solve = NULL;
172 2 : solver->iterate = redblack_omp_iterate;
173 2 : solver->apply_bc = NULL;
174 :
175 2 : return solver;
176 : }
177 : #endif
|
