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 <stdbool.h>
13 : #include <limits.h>
14 :
15 : #ifdef __cplusplus
16 : extern "C" {
17 : #endif
18 :
19 : /* ============================================================================
20 : * FACTORY FUNCTIONS
21 : *
22 : * All SIMD backends use runtime CPU detection (AVX2/NEON) via the SIMD
23 : * dispatcher. See simd/linear_solver_simd_dispatch.c for details.
24 : * ============================================================================ */
25 :
26 : /* Jacobi solvers */
27 : poisson_solver_t* create_jacobi_scalar_solver(void);
28 : poisson_solver_t* create_jacobi_simd_solver(void);
29 :
30 : /* SOR solvers */
31 : poisson_solver_t* create_sor_scalar_solver(void);
32 :
33 : /* Red-Black SOR solvers */
34 : poisson_solver_t* create_redblack_scalar_solver(void);
35 : poisson_solver_t* create_redblack_simd_solver(void);
36 :
37 : #ifdef CFD_ENABLE_OPENMP
38 : poisson_solver_t* create_redblack_omp_solver(void);
39 : poisson_solver_t* create_cg_omp_solver(void);
40 : #endif
41 :
42 : /* Conjugate Gradient solvers */
43 : poisson_solver_t* create_cg_scalar_solver(void);
44 : poisson_solver_t* create_cg_simd_solver(void);
45 :
46 : /* BiCGSTAB solvers (for non-symmetric systems) */
47 : poisson_solver_t* create_bicgstab_scalar_solver(void);
48 : poisson_solver_t* create_bicgstab_simd_solver(void);
49 :
50 : /* ============================================================================
51 : * CG ALGORITHM CONSTANTS
52 : * ============================================================================ */
53 :
54 : /**
55 : * Threshold for detecting CG breakdown (division by near-zero).
56 : * If (p, Ap) or (r, r) falls below this, the algorithm has stagnated
57 : * or encountered a singular/near-singular system.
58 : */
59 : #define CG_BREAKDOWN_THRESHOLD 1e-30
60 :
61 : /**
62 : * Macro for CG breakdown check with early return.
63 : * Used when a denominator (p_dot_Ap or r_dot_r) becomes too small.
64 : *
65 : * @param value The value to check against breakdown threshold
66 : * @param stats Pointer to stats structure (may be NULL)
67 : * @param iter Current iteration index
68 : * @param res_norm Current residual norm
69 : * @param start_time Start time for elapsed time calculation
70 : */
71 : #define CG_CHECK_BREAKDOWN(value, stats, iter, res_norm, start_time) \
72 : do { \
73 : if (fabs(value) < CG_BREAKDOWN_THRESHOLD) { \
74 : if (stats) { \
75 : (stats)->status = POISSON_STAGNATED; \
76 : (stats)->iterations = (iter) + 1; \
77 : (stats)->final_residual = (res_norm); \
78 : (stats)->elapsed_time_ms = poisson_solver_get_time_ms() - (start_time); \
79 : } \
80 : return CFD_ERROR_MAX_ITER; \
81 : } \
82 : } while (0)
83 :
84 : /* ============================================================================
85 : * BICGSTAB ALGORITHM CONSTANTS
86 : * ============================================================================ */
87 :
88 : /**
89 : * Threshold for detecting BiCGSTAB breakdown (division by near-zero).
90 : * If rho, (r_hat,v), or (t,t) falls below this, the algorithm has stagnated.
91 : */
92 : #define BICGSTAB_BREAKDOWN_THRESHOLD 1e-30
93 :
94 : /**
95 : * Convert size_t to int for OpenMP loop bounds.
96 : * OpenMP requires int loop variables, but grid dimensions are size_t.
97 : *
98 : * @param val The size_t value to convert
99 : * @return int value, or 0 on overflow (error set)
100 : */
101 : static inline int bicgstab_size_to_int(size_t val) {
102 : if (val > (size_t)INT_MAX) {
103 : cfd_set_error(CFD_ERROR_LIMIT_EXCEEDED, "Grid size exceeds INT_MAX for OpenMP loop");
104 : return 0;
105 : }
106 : return (int)val;
107 : }
108 :
109 : /* ============================================================================
110 : * SIMD BACKEND AVAILABILITY (Runtime detection)
111 : * ============================================================================ */
112 :
113 : /**
114 : * Check if SIMD backend is available at runtime.
115 : * Uses cfd_detect_simd_arch() from cpu_features.h.
116 : */
117 : bool poisson_solver_simd_backend_available(void);
118 :
119 : /**
120 : * Get the name of the detected SIMD architecture.
121 : * Returns "avx2", "neon", or "none".
122 : */
123 : const char* poisson_solver_simd_get_arch_name(void);
124 :
125 : /* ============================================================================
126 : * 3D LOOP BOUNDS HELPERS
127 : *
128 : * Centralizes the nz-dependent logic so each solver's init doesn't repeat it.
129 : * When nz==1 (2D): stride_z=0, k_start=0, k_end=1 → single k-iteration,
130 : * z-stencil terms vanish naturally.
131 : * ============================================================================ */
132 :
133 : /**
134 : * Compute 3D loop bounds from solver dimensions.
135 : */
136 183708 : static inline void poisson_solver_compute_3d_bounds(
137 : size_t nz, size_t nx, size_t ny,
138 : size_t* stride_z, size_t* k_start, size_t* k_end)
139 : {
140 183708 : *stride_z = (nz > 1) ? (nx * ny) : 0;
141 183708 : *k_start = (nz > 1) ? 1 : 0;
142 183708 : *k_end = (nz > 1) ? (nz - 1) : 1;
143 : }
144 :
145 : /**
146 : * Compute inv_dz2 safely (0.0 when dz==0, avoiding division by zero).
147 : */
148 183708 : static inline double poisson_solver_compute_inv_dz2(double dz) {
149 183708 : return (dz > 0.0) ? (1.0 / (dz * dz)) : 0.0;
150 : }
151 :
152 : /* ============================================================================
153 : * INTERNAL HELPER FUNCTIONS
154 : * ============================================================================ */
155 :
156 : /**
157 : * Common solve loop used by all iterative solvers
158 : *
159 : * Implements the iteration control, convergence checking, and statistics.
160 : *
161 : * @param solver Initialized solver
162 : * @param x Solution vector
163 : * @param x_temp Temporary buffer
164 : * @param rhs Right-hand side
165 : * @param stats Output statistics
166 : * @return CFD_SUCCESS if converged
167 : */
168 : cfd_status_t poisson_solver_solve_common(
169 : poisson_solver_t* solver,
170 : double* x,
171 : double* x_temp,
172 : const double* rhs,
173 : poisson_solver_stats_t* stats);
174 :
175 : /**
176 : * Get current time in milliseconds (platform-independent)
177 : */
178 : double poisson_solver_get_time_ms(void);
179 :
180 : #ifdef __cplusplus
181 : }
182 : #endif
183 :
184 : #endif /* CFD_LINEAR_SOLVER_INTERNAL_H */
|
