Line data Source code
1 : /**
2 : * @file linear_solver_internal.h
3 : * @brief Internal declarations for linear solver implementations
4 : *
5 : * This header is not part of the public API.
6 : */
7 :
8 : #ifndef CFD_LINEAR_SOLVER_INTERNAL_H
9 : #define CFD_LINEAR_SOLVER_INTERNAL_H
10 :
11 : #include "cfd/solvers/poisson_solver.h"
12 : #include <math.h>
13 : #include <stdbool.h>
14 : #include <limits.h>
15 :
16 : #ifndef M_PI
17 : #define M_PI 3.14159265358979323846
18 : #endif
19 :
20 : #ifdef __cplusplus
21 : extern "C" {
22 : #endif
23 :
24 : /* ============================================================================
25 : * FACTORY FUNCTIONS
26 : *
27 : * All SIMD backends use runtime CPU detection (AVX2/NEON) via the SIMD
28 : * dispatcher. See simd/linear_solver_simd_dispatch.c for details.
29 : * ============================================================================ */
30 :
31 : /* Jacobi solvers */
32 : poisson_solver_t* create_jacobi_scalar_solver(void);
33 : poisson_solver_t* create_jacobi_simd_solver(void);
34 :
35 : /* SOR solvers */
36 : poisson_solver_t* create_sor_scalar_solver(void);
37 : poisson_solver_t* create_sor_simd_solver(void);
38 :
39 : /* Red-Black SOR solvers */
40 : poisson_solver_t* create_redblack_scalar_solver(void);
41 : poisson_solver_t* create_redblack_simd_solver(void);
42 :
43 : #ifdef CFD_ENABLE_OPENMP
44 : poisson_solver_t* create_redblack_omp_solver(void);
45 : poisson_solver_t* create_cg_omp_solver(void);
46 : #endif
47 :
48 : /* Conjugate Gradient solvers */
49 : poisson_solver_t* create_cg_scalar_solver(void);
50 : poisson_solver_t* create_cg_simd_solver(void);
51 :
52 : /* BiCGSTAB solvers (for non-symmetric systems) */
53 : poisson_solver_t* create_bicgstab_scalar_solver(void);
54 : poisson_solver_t* create_bicgstab_simd_solver(void);
55 :
56 : /* ============================================================================
57 : * CG ALGORITHM CONSTANTS
58 : * ============================================================================ */
59 :
60 : /**
61 : * Threshold for detecting CG breakdown (division by near-zero).
62 : * If (p, Ap) or (r, r) falls below this, the algorithm has stagnated
63 : * or encountered a singular/near-singular system.
64 : */
65 : #define CG_BREAKDOWN_THRESHOLD 1e-30
66 :
67 : /**
68 : * Macro for CG breakdown check with early return.
69 : * Used when a denominator (p_dot_Ap or r_dot_r) becomes too small.
70 : *
71 : * @param value The value to check against breakdown threshold
72 : * @param stats Pointer to stats structure (may be NULL)
73 : * @param iter Current iteration index
74 : * @param res_norm Current residual norm
75 : * @param start_time Start time for elapsed time calculation
76 : */
77 : #define CG_CHECK_BREAKDOWN(value, stats, iter, res_norm, start_time) \
78 : do { \
79 : if (fabs(value) < CG_BREAKDOWN_THRESHOLD) { \
80 : if (stats) { \
81 : (stats)->status = POISSON_STAGNATED; \
82 : (stats)->iterations = (iter) + 1; \
83 : (stats)->final_residual = (res_norm); \
84 : (stats)->elapsed_time_ms = poisson_solver_get_time_ms() - (start_time); \
85 : } \
86 : return CFD_ERROR_MAX_ITER; \
87 : } \
88 : } while (0)
89 :
90 : /* ============================================================================
91 : * BICGSTAB ALGORITHM CONSTANTS
92 : * ============================================================================ */
93 :
94 : /**
95 : * Threshold for detecting BiCGSTAB breakdown (division by near-zero).
96 : * If rho, (r_hat,v), or (t,t) falls below this, the algorithm has stagnated.
97 : */
98 : #define BICGSTAB_BREAKDOWN_THRESHOLD 1e-30
99 :
100 : /**
101 : * Convert size_t to int for OpenMP loop bounds.
102 : * OpenMP requires int loop variables, but grid dimensions are size_t.
103 : *
104 : * @param val The size_t value to convert
105 : * @return int value, or 0 on overflow (error set)
106 : */
107 : static inline int bicgstab_size_to_int(size_t val) {
108 : if (val > (size_t)INT_MAX) {
109 : cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop");
110 : return 0;
111 : }
112 : return (int)val;
113 : }
114 :
115 : /* ============================================================================
116 : * SIMD BACKEND AVAILABILITY (Runtime detection)
117 : * ============================================================================ */
118 :
119 : /**
120 : * Check if SIMD backend is available at runtime.
121 : * Uses cfd_detect_simd_arch() from cpu_features.h.
122 : */
123 : bool poisson_solver_simd_backend_available(void);
124 :
125 : /**
126 : * Get the name of the detected SIMD architecture.
127 : * Returns "avx2", "neon", or "none".
128 : */
129 : const char* poisson_solver_simd_get_arch_name(void);
130 :
131 : /* ============================================================================
132 : * 3D LOOP BOUNDS HELPERS
133 : *
134 : * Centralizes the nz-dependent logic so each solver's init doesn't repeat it.
135 : * When nz==1 (2D): stride_z=0, k_start=0, k_end=1 → single k-iteration,
136 : * z-stencil terms vanish naturally.
137 : * ============================================================================ */
138 :
139 : /**
140 : * Compute 3D loop bounds from solver dimensions.
141 : */
142 185804 : static inline void poisson_solver_compute_3d_bounds(
143 : size_t nz, size_t nx, size_t ny,
144 : size_t* stride_z, size_t* k_start, size_t* k_end)
145 : {
146 185804 : *stride_z = (nz > 1) ? (nx * ny) : 0;
147 185804 : *k_start = (nz > 1) ? 1 : 0;
148 185804 : *k_end = (nz > 1) ? (nz - 1) : 1;
149 : }
150 :
151 : /**
152 : * Compute inv_dz2 safely (0.0 when dz==0, avoiding division by zero).
153 : */
154 185816 : static inline double poisson_solver_compute_inv_dz2(double dz) {
155 185804 : return (dz > 0.0) ? (1.0 / (dz * dz)) : 0.0;
156 : }
157 :
158 : /* ============================================================================
159 : * OPTIMAL SOR OMEGA COMPUTATION
160 : * ============================================================================ */
161 :
162 : /**
163 : * Compute optimal SOR relaxation parameter for a 2D/3D Laplacian.
164 : *
165 : * Uses the Jacobi spectral radius for a rectangular grid:
166 : * 2D: rho_J = [cos(pi/(nx-1))/dx^2 + cos(pi/(ny-1))/dy^2]
167 : * / [1/dx^2 + 1/dy^2]
168 : * 3D: rho_J = [cos(pi/(nx-1))/dx^2 + cos(pi/(ny-1))/dy^2
169 : * + cos(pi/(nz-1))/dz^2]
170 : * / [1/dx^2 + 1/dy^2 + 1/dz^2]
171 : * omega_opt = 2 / (1 + sqrt(1 - rho_J^2))
172 : *
173 : * When nz <= 1 or dz <= 0, the z-component is ignored and the 2D formula
174 : * is used.
175 : */
176 12 : static inline double poisson_solver_compute_optimal_omega(
177 : size_t nx, size_t ny, size_t nz,
178 : double dx, double dy, double dz)
179 : {
180 12 : double inv_dx2 = 1.0 / (dx * dx);
181 12 : double inv_dy2 = 1.0 / (dy * dy);
182 12 : double inv_dz2 = poisson_solver_compute_inv_dz2(dz);
183 :
184 12 : double num = cos(M_PI / (double)(nx - 1)) * inv_dx2
185 12 : + cos(M_PI / (double)(ny - 1)) * inv_dy2;
186 12 : double denom = inv_dx2 + inv_dy2;
187 :
188 12 : if (nz > 1 && inv_dz2 > 0.0) {
189 0 : num += cos(M_PI / (double)(nz - 1)) * inv_dz2;
190 0 : denom += inv_dz2;
191 : }
192 :
193 12 : double rho_j = num / denom;
194 12 : return 2.0 / (1.0 + sqrt(1.0 - (rho_j * rho_j)));
195 : }
196 :
197 : /**
198 : * Resolve omega: if omega <= 0.0 (auto sentinel), compute optimal value.
199 : * Otherwise return the user-specified omega as-is.
200 : *
201 : * Supports both 2D (nz <= 1 or dz <= 0) and 3D (nz > 1, dz > 0) problems.
202 : */
203 46 : static inline double poisson_solver_resolve_omega(
204 : double omega,
205 : size_t nx, size_t ny, size_t nz,
206 : double dx, double dy, double dz)
207 : {
208 46 : if (omega <= 0.0) {
209 12 : return poisson_solver_compute_optimal_omega(nx, ny, nz, dx, dy, dz);
210 : }
211 : return omega;
212 : }
213 :
214 : /* ============================================================================
215 : * INTERNAL HELPER FUNCTIONS
216 : * ============================================================================ */
217 :
218 : /**
219 : * Common solve loop used by all iterative solvers
220 : *
221 : * Implements the iteration control, convergence checking, and statistics.
222 : *
223 : * @param solver Initialized solver
224 : * @param x Solution vector
225 : * @param x_temp Temporary buffer
226 : * @param rhs Right-hand side
227 : * @param stats Output statistics
228 : * @return CFD_SUCCESS if converged
229 : */
230 : cfd_status_t poisson_solver_solve_common(
231 : poisson_solver_t* solver,
232 : double* x,
233 : double* x_temp,
234 : const double* rhs,
235 : poisson_solver_stats_t* stats);
236 :
237 : /**
238 : * Get current time in milliseconds (platform-independent)
239 : */
240 : double poisson_solver_get_time_ms(void);
241 :
242 : #ifdef __cplusplus
243 : }
244 : #endif
245 :
246 : #endif /* CFD_LINEAR_SOLVER_INTERNAL_H */
|
