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 3 : 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 3 : redblack_omp_context_t* ctx = (redblack_omp_context_t*)cfd_calloc(1, sizeof(redblack_omp_context_t));
52 3 : if (!ctx) {
53 : return CFD_ERROR_NOMEM;
54 : }
55 :
56 3 : ctx->dx2 = dx * dx;
57 3 : ctx->dy2 = dy * dy;
58 3 : ctx->inv_dz2 = poisson_solver_compute_inv_dz2(dz);
59 3 : poisson_solver_compute_3d_bounds(nz, nx, ny, &ctx->stride_z, &ctx->k_start, &ctx->k_end);
60 3 : double factor = 2.0 * (1.0 / ctx->dx2 + 1.0 / ctx->dy2 + ctx->inv_dz2);
61 3 : ctx->inv_factor = 1.0 / factor;
62 3 : ctx->omega = poisson_solver_resolve_omega(
63 : params ? params->omega : 0.0, nx, ny, nz, dx, dy, dz);
64 3 : ctx->initialized = 1;
65 :
66 3 : solver->context = ctx;
67 3 : return CFD_SUCCESS;
68 : }
69 :
70 5 : static void redblack_omp_destroy(poisson_solver_t* solver) {
71 5 : if (solver && solver->context) {
72 3 : cfd_free(solver->context);
73 3 : solver->context = NULL;
74 : }
75 5 : }
76 :
77 327 : static cfd_status_t redblack_omp_iterate(
78 : poisson_solver_t* solver,
79 : double* x,
80 : double* x_temp,
81 : const double* rhs,
82 : double* residual)
83 : {
84 327 : (void)x_temp;
85 :
86 327 : redblack_omp_context_t* ctx = (redblack_omp_context_t*)solver->context;
87 327 : size_t nx = solver->nx;
88 327 : size_t ny = solver->ny;
89 327 : double dx2 = ctx->dx2;
90 327 : double dy2 = ctx->dy2;
91 327 : double inv_dz2 = ctx->inv_dz2;
92 327 : double inv_factor = ctx->inv_factor;
93 327 : double omega = ctx->omega;
94 327 : size_t stride_z = ctx->stride_z;
95 327 : size_t k_start = ctx->k_start;
96 327 : size_t k_end = ctx->k_end;
97 :
98 : /* Red sweep (parallel) */
99 2852 : for (size_t k = k_start; k < k_end; k++) {
100 2525 : int j;
101 2525 : #pragma omp parallel for schedule(static)
102 : for (j = 1; j < (int)ny - 1; j++) {
103 : size_t i_start = ((j + k) % 2 == 0) ? 1 : 2;
104 : size_t i;
105 : for (i = i_start; i < nx - 1; i += 2) {
106 : size_t idx = k * stride_z + IDX_2D(i, (size_t)j, nx);
107 :
108 : double p_new = -(rhs[idx]
109 : - (x[idx + 1] + x[idx - 1]) / dx2
110 : - (x[idx + nx] + x[idx - nx]) / dy2
111 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
112 : ) * inv_factor;
113 :
114 : x[idx] = x[idx] + omega * (p_new - x[idx]);
115 : }
116 : }
117 : }
118 :
119 : /* Black sweep (parallel) */
120 2852 : for (size_t k = k_start; k < k_end; k++) {
121 2525 : int j;
122 2525 : #pragma omp parallel for schedule(static)
123 : for (j = 1; j < (int)ny - 1; j++) {
124 : size_t i_start = ((j + k) % 2 == 0) ? 2 : 1;
125 : size_t i;
126 : for (i = i_start; i < nx - 1; i += 2) {
127 : size_t idx = k * stride_z + IDX_2D(i, (size_t)j, nx);
128 :
129 : double p_new = -(rhs[idx]
130 : - (x[idx + 1] + x[idx - 1]) / dx2
131 : - (x[idx + nx] + x[idx - nx]) / dy2
132 : - (x[idx + stride_z] + x[idx - stride_z]) * inv_dz2
133 : ) * inv_factor;
134 :
135 : x[idx] = x[idx] + omega * (p_new - x[idx]);
136 : }
137 : }
138 : }
139 :
140 : /* Apply boundary conditions */
141 327 : poisson_solver_apply_bc(solver, x);
142 :
143 : /* Compute residual if requested */
144 327 : if (residual) {
145 327 : *residual = poisson_solver_compute_residual(solver, x, rhs);
146 : }
147 :
148 327 : return CFD_SUCCESS;
149 : }
150 :
151 : #endif /* CFD_ENABLE_OPENMP */
152 :
153 : /* ============================================================================
154 : * FACTORY FUNCTION
155 : * ============================================================================ */
156 :
157 : #ifdef CFD_ENABLE_OPENMP
158 5 : poisson_solver_t* create_redblack_omp_solver(void) {
159 5 : poisson_solver_t* solver = (poisson_solver_t*)cfd_calloc(1, sizeof(poisson_solver_t));
160 5 : if (!solver) {
161 : return NULL;
162 : }
163 :
164 5 : solver->name = POISSON_SOLVER_TYPE_REDBLACK_OMP;
165 5 : solver->description = "Red-Black SOR iteration (OpenMP parallel)";
166 5 : solver->method = POISSON_METHOD_REDBLACK_SOR;
167 5 : solver->backend = POISSON_BACKEND_OMP;
168 5 : solver->params = poisson_solver_params_default();
169 :
170 5 : solver->init = redblack_omp_init;
171 5 : solver->destroy = redblack_omp_destroy;
172 5 : solver->solve = NULL;
173 5 : solver->iterate = redblack_omp_iterate;
174 5 : solver->apply_bc = NULL;
175 :
176 5 : return solver;
177 : }
178 : #endif
|
