Line data Source code
1 : /**
2 : * @file solver_rk2_avx2.c
3 : * @brief RK2 (Heun's method) - AVX2 + OpenMP implementation
4 : *
5 : * AVX2-only backend. Returns CFD_ERROR_UNSUPPORTED if AVX2 is unavailable.
6 : * Uses persistent aligned buffers in context to avoid per-step allocation.
7 : *
8 : * SIMD strategy per row (fixed j, k):
9 : * i = 1 : scalar (periodic il wrap)
10 : * i = 2..nx-4 : AVX2 vectorized (interior, no wrapping)
11 : * i = nx-2 : scalar (periodic ir wrap)
12 : * Remainder between AVX2 end and nx-2: scalar (no wrapping)
13 : *
14 : * Algorithm: identical to scalar/OMP RK2.
15 : * k1 = RHS(Q^n)
16 : * Q_pred = Q^n + dt*k1
17 : * k2 = RHS(Q_pred)
18 : * Q^{n+1} = Q^n + (dt/2)*(k1 + k2)
19 : * BCs applied to final state only.
20 : *
21 : * Branch-free 3D: when nz==1, stride_z=0 and inv_2dz/inv_dz2=0.0 cause all
22 : * z-terms to vanish, producing bit-identical results to the 2D code path.
23 : */
24 :
25 : #define _POSIX_C_SOURCE 200809L
26 : #define _ISOC11_SOURCE
27 : #define _USE_MATH_DEFINES
28 :
29 : #include "cfd/core/cfd_status.h"
30 : #include "cfd/core/grid.h"
31 : #include "cfd/core/indexing.h"
32 : #include "cfd/core/logging.h"
33 : #include "cfd/core/memory.h"
34 : #include "cfd/solvers/navier_stokes_solver.h"
35 : #include "cfd/solvers/energy_solver.h"
36 : #include "../../energy/energy_solver_internal.h"
37 :
38 : #include <math.h>
39 : #include <stddef.h>
40 : #include <stdio.h>
41 : #include <string.h>
42 :
43 : #ifndef M_PI
44 : #define M_PI 3.14159265358979323846
45 : #endif
46 :
47 : #ifdef _OPENMP
48 : #include <omp.h>
49 : #endif
50 :
51 : #if defined(CFD_HAS_AVX2)
52 : #include <immintrin.h>
53 : #define USE_AVX2 1
54 : #else
55 : #define USE_AVX2 0
56 : #endif
57 :
58 : /* Physical stability limits (identical to scalar/OMP RK2) */
59 : #define MAX_DERIVATIVE_LIMIT 100.0
60 : #define MAX_SECOND_DERIVATIVE_LIMIT 1000.0
61 : #define MAX_VELOCITY_LIMIT 100.0
62 : #define MAX_DIVERGENCE_LIMIT 10.0
63 : #define PRESSURE_UPDATE_FACTOR 0.1
64 :
65 : /* ============================================================================
66 : * CONTEXT
67 : * ============================================================================ */
68 :
69 : typedef struct {
70 : double* k1_u; double* k1_v; double* k1_w; double* k1_p;
71 : double* k2_u; double* k2_v; double* k2_w; double* k2_p;
72 : double* u0; double* v0; double* w0; double* p0;
73 : double* T_ws; /* Reusable scratch for the energy step (avoids per-step alloc) */
74 : double* dx_inv; /* 1/(2*dx[i]) for i = 0..nx-1 */
75 : double* dy_inv; /* 1/(2*dy[j]) for j = 0..ny-1 */
76 : size_t nx, ny, nz;
77 : size_t stride_z;
78 : size_t k_start, k_end;
79 : double inv_2dz, inv_dz2;
80 : int initialized;
81 : } rk2_avx2_context_t;
82 :
83 : /* ============================================================================
84 : * PUBLIC API DECLARATIONS
85 : * ============================================================================ */
86 :
87 : cfd_status_t rk2_avx2_init(ns_solver_t* solver, const grid* g,
88 : const ns_solver_params_t* params);
89 : void rk2_avx2_destroy(ns_solver_t* solver);
90 : cfd_status_t rk2_avx2_step(ns_solver_t* solver, flow_field* field, const grid* g,
91 : const ns_solver_params_t* params, ns_solver_stats_t* stats);
92 :
93 : /* ============================================================================
94 : * SCALAR POINT HELPER
95 : * Handles a single interior grid point given explicit stencil indices.
96 : * Used for periodic-edge columns (i=1 with il wrap, i=nx-2 with ir wrap)
97 : * and for scalar remainder after the AVX2 group.
98 : * ============================================================================ */
99 :
100 : static void rk2_rhs_point(
101 : const double* u, const double* v, const double* w,
102 : const double* p, const double* rho, const double* T,
103 : double* rhs_u, double* rhs_v, double* rhs_w, double* rhs_p,
104 : const grid* g, const ns_solver_params_t* params,
105 : size_t idx, size_t il, size_t ir, size_t jd, size_t ju,
106 : size_t kd, size_t ku,
107 : double inv_2dz, double inv_dz2,
108 : size_t i, size_t j, size_t k, size_t nz,
109 : int iter, double dt)
110 : {
111 : if (rho[idx] <= 1e-10 || fabs(g->dx[i]) < 1e-10 || fabs(g->dy[j]) < 1e-10) {
112 : rhs_u[idx] = 0.0;
113 : rhs_v[idx] = 0.0;
114 : rhs_w[idx] = 0.0;
115 : rhs_p[idx] = 0.0;
116 : return;
117 : }
118 :
119 : double dx = g->dx[i], dy = g->dy[j];
120 : double dx2 = dx * dx, dy2 = dy * dy;
121 :
122 : double du_dx = (u[ir] - u[il]) / (2.0 * dx);
123 : double du_dy = (u[ju] - u[jd]) / (2.0 * dy);
124 : double du_dz = (u[ku] - u[kd]) * inv_2dz;
125 : double dv_dx = (v[ir] - v[il]) / (2.0 * dx);
126 : double dv_dy = (v[ju] - v[jd]) / (2.0 * dy);
127 : double dv_dz = (v[ku] - v[kd]) * inv_2dz;
128 : double dw_dx = (w[ir] - w[il]) / (2.0 * dx);
129 : double dw_dy = (w[ju] - w[jd]) / (2.0 * dy);
130 : double dw_dz = (w[ku] - w[kd]) * inv_2dz;
131 : double dp_dx = (p[ir] - p[il]) / (2.0 * dx);
132 : double dp_dy = (p[ju] - p[jd]) / (2.0 * dy);
133 : double dp_dz = (p[ku] - p[kd]) * inv_2dz;
134 :
135 : double d2u_dx2 = (u[ir] - 2.0 * u[idx] + u[il]) / dx2;
136 : double d2u_dy2 = (u[ju] - 2.0 * u[idx] + u[jd]) / dy2;
137 : double d2u_dz2 = (u[ku] - 2.0 * u[idx] + u[kd]) * inv_dz2;
138 : double d2v_dx2 = (v[ir] - 2.0 * v[idx] + v[il]) / dx2;
139 : double d2v_dy2 = (v[ju] - 2.0 * v[idx] + v[jd]) / dy2;
140 : double d2v_dz2 = (v[ku] - 2.0 * v[idx] + v[kd]) * inv_dz2;
141 : double d2w_dx2 = (w[ir] - 2.0 * w[idx] + w[il]) / dx2;
142 : double d2w_dy2 = (w[ju] - 2.0 * w[idx] + w[jd]) / dy2;
143 : double d2w_dz2 = (w[ku] - 2.0 * w[idx] + w[kd]) * inv_dz2;
144 :
145 : double nu = params->mu / fmax(rho[idx], 1e-10);
146 : nu = fmin(nu, 1.0);
147 :
148 : du_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dx));
149 : du_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dy));
150 : du_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dz));
151 : dv_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dx));
152 : dv_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dy));
153 : dv_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dz));
154 : dw_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dx));
155 : dw_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dy));
156 : dw_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dz));
157 : dp_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dx));
158 : dp_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dy));
159 : dp_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dz));
160 : d2u_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dx2));
161 : d2u_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dy2));
162 : d2u_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dz2));
163 : d2v_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dx2));
164 : d2v_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dy2));
165 : d2v_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dz2));
166 : d2w_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dx2));
167 : d2w_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dy2));
168 : d2w_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dz2));
169 :
170 : double source_u = 0.0, source_v = 0.0, source_w = 0.0;
171 : double z_coord = (nz > 1 && g->z) ? g->z[k] : 0.0;
172 : compute_source_terms(g->x[i], g->y[j], z_coord, iter, dt, params,
173 : &source_u, &source_v, &source_w);
174 :
175 : /* Boussinesq buoyancy source (no-op when beta == 0) */
176 : if (T) {
177 : energy_compute_buoyancy(T[idx], params, &source_u, &source_v, &source_w);
178 : }
179 :
180 : rhs_u[idx] = -(u[idx] * du_dx) - (v[idx] * du_dy) - (w[idx] * du_dz)
181 : - dp_dx / rho[idx]
182 : + (nu * (d2u_dx2 + d2u_dy2 + d2u_dz2))
183 : + source_u;
184 :
185 : rhs_v[idx] = -(u[idx] * dv_dx) - (v[idx] * dv_dy) - (w[idx] * dv_dz)
186 : - dp_dy / rho[idx]
187 : + (nu * (d2v_dx2 + d2v_dy2 + d2v_dz2))
188 : + source_v;
189 :
190 : rhs_w[idx] = -(u[idx] * dw_dx) - (v[idx] * dw_dy) - (w[idx] * dw_dz)
191 : - dp_dz / rho[idx]
192 : + (nu * (d2w_dx2 + d2w_dy2 + d2w_dz2))
193 : + source_w;
194 :
195 : double divergence = du_dx + dv_dy + dw_dz;
196 : divergence = fmax(-MAX_DIVERGENCE_LIMIT, fmin(MAX_DIVERGENCE_LIMIT, divergence));
197 : rhs_p[idx] = -PRESSURE_UPDATE_FACTOR * rho[idx] * divergence;
198 : }
199 :
200 : /* ============================================================================
201 : * AVX2 ROW-LEVEL RHS (compiled only when AVX2 is available)
202 : * ============================================================================ */
203 :
204 : #if USE_AVX2
205 :
206 : static inline __m256d avx2_clamp(__m256d x, __m256d lo, __m256d hi) {
207 : return _mm256_max_pd(lo, _mm256_min_pd(hi, x));
208 : }
209 :
210 : /*
211 : * Compute RHS for one row j at z-plane k, dispatching:
212 : * i=1 : scalar (periodic il wrap)
213 : * i=2..nx-4 : AVX2 (4-wide, no wrapping)
214 : * i=nx-2 : scalar (periodic ir wrap), plus any scalar remainder
215 : */
216 : static void compute_rhs_row(
217 : const double* u, const double* v, const double* w,
218 : const double* p, const double* rho, const double* T,
219 : double* rhs_u, double* rhs_v, double* rhs_w, double* rhs_p,
220 : const rk2_avx2_context_t* ctx, const grid* g,
221 : const ns_solver_params_t* params, size_t j, size_t k, int iter, double dt)
222 : {
223 : size_t nx = ctx->nx;
224 : size_t ny = ctx->ny;
225 : size_t nz = ctx->nz;
226 : size_t stride_z = ctx->stride_z;
227 : double inv_2dz = ctx->inv_2dz;
228 : double inv_dz2 = ctx->inv_dz2;
229 :
230 : size_t k_off = k * stride_z;
231 : size_t j_off = k_off + j * nx;
232 :
233 : /* Degenerate dy: zero RHS for entire row (matches scalar path) */
234 : if (fabs(g->dy[j]) < 1e-10) {
235 : return; /* rhs arrays already memset to 0 before this call */
236 : }
237 :
238 : /* j-direction stencil row offsets: uniform across entire row */
239 : size_t jd_row = (j > 1) ? k_off + (j - 1) * nx : k_off + (ny - 2) * nx;
240 : size_t ju_row = (j < ny - 2) ? k_off + (j + 1) * nx : k_off + nx;
241 :
242 : /* z-direction stencil plane offsets: uniform across entire row.
243 : * When nz==1: stride_z=0, so kd_off=ku_off=0 and z-terms vanish. */
244 : size_t kd_off = (k > 1) ? (k - 1) * stride_z : (nz - 2) * stride_z;
245 : size_t ku_off = (k < nz - 2) ? (k + 1) * stride_z : 1 * stride_z;
246 :
247 : /* --- Scalar: i = 1 (il wraps periodically) --- */
248 : {
249 : size_t si = 1;
250 : size_t idx = j_off + si;
251 : size_t il = j_off + (nx - 2); /* periodic wrap */
252 : size_t ir = idx + 1;
253 : size_t kd = kd_off + j * nx + si;
254 : size_t ku = ku_off + j * nx + si;
255 : rk2_rhs_point(u, v, w, p, rho, T, rhs_u, rhs_v, rhs_w, rhs_p, g, params,
256 : idx, il, ir, jd_row + si, ju_row + si, kd, ku,
257 : inv_2dz, inv_dz2, si, j, k, nz, iter, dt);
258 : }
259 :
260 : /* --- AVX2: i = 2 while i+3 <= nx-3 (interior, no wrapping needed) --- */
261 : {
262 : __m256d dy_inv_v = _mm256_set1_pd(ctx->dy_inv[j]);
263 : /* 1/dy^2 = 4 * (1/(2*dy))^2, derived from pre-guarded dy_inv */
264 : __m256d dy2_inv_v = _mm256_set1_pd(4.0 * ctx->dy_inv[j] * ctx->dy_inv[j]);
265 : __m256d dz_inv_v = _mm256_set1_pd(inv_2dz);
266 : __m256d dz2_inv_v = _mm256_set1_pd(inv_dz2);
267 : __m256d two = _mm256_set1_pd(2.0);
268 : __m256d four = _mm256_set1_pd(4.0);
269 : __m256d max_d1 = _mm256_set1_pd( MAX_DERIVATIVE_LIMIT);
270 : __m256d min_d1 = _mm256_set1_pd(-MAX_DERIVATIVE_LIMIT);
271 : __m256d max_d2 = _mm256_set1_pd( MAX_SECOND_DERIVATIVE_LIMIT);
272 : __m256d min_d2 = _mm256_set1_pd(-MAX_SECOND_DERIVATIVE_LIMIT);
273 : __m256d max_div = _mm256_set1_pd( MAX_DIVERGENCE_LIMIT);
274 : __m256d min_div = _mm256_set1_pd(-MAX_DIVERGENCE_LIMIT);
275 : __m256d mu_v = _mm256_set1_pd(params->mu);
276 : __m256d eps = _mm256_set1_pd(1e-10);
277 : __m256d one = _mm256_set1_pd(1.0);
278 : __m256d neg_puf = _mm256_set1_pd(-PRESSURE_UPDATE_FACTOR);
279 :
280 : ptrdiff_t i;
281 : for (i = 2; i + 3 <= (ptrdiff_t)(nx - 3); i += 4) {
282 : size_t idx = j_off + (size_t)i;
283 :
284 : /* Load u stencil */
285 : __m256d u_c = _mm256_loadu_pd(&u[idx]);
286 : __m256d u_il = _mm256_loadu_pd(&u[idx - 1]);
287 : __m256d u_ir = _mm256_loadu_pd(&u[idx + 1]);
288 : __m256d u_jd = _mm256_loadu_pd(&u[jd_row + (size_t)i]);
289 : __m256d u_ju = _mm256_loadu_pd(&u[ju_row + (size_t)i]);
290 : __m256d u_kd = _mm256_loadu_pd(&u[kd_off + j * nx + (size_t)i]);
291 : __m256d u_ku = _mm256_loadu_pd(&u[ku_off + j * nx + (size_t)i]);
292 :
293 : /* Load v stencil */
294 : __m256d v_c = _mm256_loadu_pd(&v[idx]);
295 : __m256d v_il = _mm256_loadu_pd(&v[idx - 1]);
296 : __m256d v_ir = _mm256_loadu_pd(&v[idx + 1]);
297 : __m256d v_jd = _mm256_loadu_pd(&v[jd_row + (size_t)i]);
298 : __m256d v_ju = _mm256_loadu_pd(&v[ju_row + (size_t)i]);
299 : __m256d v_kd = _mm256_loadu_pd(&v[kd_off + j * nx + (size_t)i]);
300 : __m256d v_ku = _mm256_loadu_pd(&v[ku_off + j * nx + (size_t)i]);
301 :
302 : /* Load w stencil */
303 : __m256d w_c = _mm256_loadu_pd(&w[idx]);
304 : __m256d w_il = _mm256_loadu_pd(&w[idx - 1]);
305 : __m256d w_ir = _mm256_loadu_pd(&w[idx + 1]);
306 : __m256d w_jd = _mm256_loadu_pd(&w[jd_row + (size_t)i]);
307 : __m256d w_ju = _mm256_loadu_pd(&w[ju_row + (size_t)i]);
308 : __m256d w_kd = _mm256_loadu_pd(&w[kd_off + j * nx + (size_t)i]);
309 : __m256d w_ku = _mm256_loadu_pd(&w[ku_off + j * nx + (size_t)i]);
310 :
311 : /* Load p stencil */
312 : __m256d p_il = _mm256_loadu_pd(&p[idx - 1]);
313 : __m256d p_ir = _mm256_loadu_pd(&p[idx + 1]);
314 : __m256d p_jd = _mm256_loadu_pd(&p[jd_row + (size_t)i]);
315 : __m256d p_ju = _mm256_loadu_pd(&p[ju_row + (size_t)i]);
316 : __m256d p_kd = _mm256_loadu_pd(&p[kd_off + j * nx + (size_t)i]);
317 : __m256d p_ku = _mm256_loadu_pd(&p[ku_off + j * nx + (size_t)i]);
318 :
319 : /* Density: clamp for safe division, but track validity */
320 : __m256d rho_raw = _mm256_loadu_pd(&rho[idx]);
321 : __m256d rho_valid = _mm256_cmp_pd(rho_raw, eps, _CMP_GT_OQ);
322 : __m256d rho_c = _mm256_max_pd(rho_raw, eps);
323 : __m256d rho_inv = _mm256_div_pd(one, rho_c);
324 :
325 : /* dx stencil: dx_inv[i] = 1/(2*dx[i]), 0 when dx < 1e-10 */
326 : __m256d dx_inv_v = _mm256_loadu_pd(&ctx->dx_inv[(size_t)i]);
327 : __m256d dx_valid = _mm256_cmp_pd(dx_inv_v, _mm256_setzero_pd(), _CMP_GT_OQ);
328 :
329 : /* Combined validity: rho > eps AND dx > 0 (dy already checked at row level) */
330 : __m256d valid = _mm256_and_pd(rho_valid, dx_valid);
331 : /* 1/dx^2 = 4 * (1/(2*dx))^2 */
332 : __m256d dx2_inv_v = _mm256_mul_pd(four, _mm256_mul_pd(dx_inv_v, dx_inv_v));
333 :
334 : /* First derivatives — u */
335 : __m256d du_dx = avx2_clamp(
336 : _mm256_mul_pd(_mm256_sub_pd(u_ir, u_il), dx_inv_v), min_d1, max_d1);
337 : __m256d du_dy = avx2_clamp(
338 : _mm256_mul_pd(_mm256_sub_pd(u_ju, u_jd), dy_inv_v), min_d1, max_d1);
339 : __m256d du_dz = avx2_clamp(
340 : _mm256_mul_pd(_mm256_sub_pd(u_ku, u_kd), dz_inv_v), min_d1, max_d1);
341 :
342 : /* First derivatives — v */
343 : __m256d dv_dx = avx2_clamp(
344 : _mm256_mul_pd(_mm256_sub_pd(v_ir, v_il), dx_inv_v), min_d1, max_d1);
345 : __m256d dv_dy = avx2_clamp(
346 : _mm256_mul_pd(_mm256_sub_pd(v_ju, v_jd), dy_inv_v), min_d1, max_d1);
347 : __m256d dv_dz = avx2_clamp(
348 : _mm256_mul_pd(_mm256_sub_pd(v_ku, v_kd), dz_inv_v), min_d1, max_d1);
349 :
350 : /* First derivatives — w */
351 : __m256d dw_dx = avx2_clamp(
352 : _mm256_mul_pd(_mm256_sub_pd(w_ir, w_il), dx_inv_v), min_d1, max_d1);
353 : __m256d dw_dy = avx2_clamp(
354 : _mm256_mul_pd(_mm256_sub_pd(w_ju, w_jd), dy_inv_v), min_d1, max_d1);
355 : __m256d dw_dz = avx2_clamp(
356 : _mm256_mul_pd(_mm256_sub_pd(w_ku, w_kd), dz_inv_v), min_d1, max_d1);
357 :
358 : /* Pressure gradients */
359 : __m256d dp_dx = avx2_clamp(
360 : _mm256_mul_pd(_mm256_sub_pd(p_ir, p_il), dx_inv_v), min_d1, max_d1);
361 : __m256d dp_dy = avx2_clamp(
362 : _mm256_mul_pd(_mm256_sub_pd(p_ju, p_jd), dy_inv_v), min_d1, max_d1);
363 : __m256d dp_dz = avx2_clamp(
364 : _mm256_mul_pd(_mm256_sub_pd(p_ku, p_kd), dz_inv_v), min_d1, max_d1);
365 :
366 : /* Second derivatives — u */
367 : __m256d d2u_dx2 = avx2_clamp(
368 : _mm256_mul_pd(
369 : _mm256_sub_pd(_mm256_add_pd(u_ir, u_il), _mm256_mul_pd(two, u_c)),
370 : dx2_inv_v),
371 : min_d2, max_d2);
372 : __m256d d2u_dy2 = avx2_clamp(
373 : _mm256_mul_pd(
374 : _mm256_sub_pd(_mm256_add_pd(u_ju, u_jd), _mm256_mul_pd(two, u_c)),
375 : dy2_inv_v),
376 : min_d2, max_d2);
377 : __m256d d2u_dz2 = avx2_clamp(
378 : _mm256_mul_pd(
379 : _mm256_sub_pd(_mm256_add_pd(u_ku, u_kd), _mm256_mul_pd(two, u_c)),
380 : dz2_inv_v),
381 : min_d2, max_d2);
382 :
383 : /* Second derivatives — v */
384 : __m256d d2v_dx2 = avx2_clamp(
385 : _mm256_mul_pd(
386 : _mm256_sub_pd(_mm256_add_pd(v_ir, v_il), _mm256_mul_pd(two, v_c)),
387 : dx2_inv_v),
388 : min_d2, max_d2);
389 : __m256d d2v_dy2 = avx2_clamp(
390 : _mm256_mul_pd(
391 : _mm256_sub_pd(_mm256_add_pd(v_ju, v_jd), _mm256_mul_pd(two, v_c)),
392 : dy2_inv_v),
393 : min_d2, max_d2);
394 : __m256d d2v_dz2 = avx2_clamp(
395 : _mm256_mul_pd(
396 : _mm256_sub_pd(_mm256_add_pd(v_ku, v_kd), _mm256_mul_pd(two, v_c)),
397 : dz2_inv_v),
398 : min_d2, max_d2);
399 :
400 : /* Second derivatives — w */
401 : __m256d d2w_dx2 = avx2_clamp(
402 : _mm256_mul_pd(
403 : _mm256_sub_pd(_mm256_add_pd(w_ir, w_il), _mm256_mul_pd(two, w_c)),
404 : dx2_inv_v),
405 : min_d2, max_d2);
406 : __m256d d2w_dy2 = avx2_clamp(
407 : _mm256_mul_pd(
408 : _mm256_sub_pd(_mm256_add_pd(w_ju, w_jd), _mm256_mul_pd(two, w_c)),
409 : dy2_inv_v),
410 : min_d2, max_d2);
411 : __m256d d2w_dz2 = avx2_clamp(
412 : _mm256_mul_pd(
413 : _mm256_sub_pd(_mm256_add_pd(w_ku, w_kd), _mm256_mul_pd(two, w_c)),
414 : dz2_inv_v),
415 : min_d2, max_d2);
416 :
417 : /* Kinematic viscosity: nu = min(mu/rho, 1.0) */
418 : __m256d nu = _mm256_min_pd(one, _mm256_mul_pd(mu_v, rho_inv));
419 :
420 : /* Source terms: computed scalar per lane and packed */
421 : double src_u_arr[4] = {0.0, 0.0, 0.0, 0.0};
422 : double src_v_arr[4] = {0.0, 0.0, 0.0, 0.0};
423 : double src_w_arr[4] = {0.0, 0.0, 0.0, 0.0};
424 : for (int lane = 0; lane < 4; lane++) {
425 : double z_coord = (nz > 1 && g->z) ? g->z[k] : 0.0;
426 : compute_source_terms(g->x[(size_t)i + (size_t)lane], g->y[j],
427 : z_coord, iter, dt, params,
428 : &src_u_arr[lane], &src_v_arr[lane],
429 : &src_w_arr[lane]);
430 : /* Boussinesq buoyancy source (no-op when beta == 0) */
431 : if (T) {
432 : energy_compute_buoyancy(T[j_off + (size_t)i + (size_t)lane], params,
433 : &src_u_arr[lane], &src_v_arr[lane],
434 : &src_w_arr[lane]);
435 : }
436 : }
437 : __m256d src_u_v = _mm256_loadu_pd(src_u_arr);
438 : __m256d src_v_v = _mm256_loadu_pd(src_v_arr);
439 : __m256d src_w_v = _mm256_loadu_pd(src_w_arr);
440 :
441 : /* RHS u: nu*(d2u_dx2+d2u_dy2+d2u_dz2) - dp_dx/rho - conv_u + src_u */
442 : __m256d conv_u = _mm256_add_pd(
443 : _mm256_add_pd(
444 : _mm256_mul_pd(u_c, du_dx),
445 : _mm256_mul_pd(v_c, du_dy)),
446 : _mm256_mul_pd(w_c, du_dz));
447 : __m256d rhs_u_v = _mm256_add_pd(
448 : _mm256_sub_pd(
449 : _mm256_sub_pd(
450 : _mm256_mul_pd(nu, _mm256_add_pd(_mm256_add_pd(d2u_dx2, d2u_dy2), d2u_dz2)),
451 : _mm256_mul_pd(dp_dx, rho_inv)),
452 : conv_u),
453 : src_u_v);
454 :
455 : /* RHS v: nu*(d2v_dx2+d2v_dy2+d2v_dz2) - dp_dy/rho - conv_v + src_v */
456 : __m256d conv_v = _mm256_add_pd(
457 : _mm256_add_pd(
458 : _mm256_mul_pd(u_c, dv_dx),
459 : _mm256_mul_pd(v_c, dv_dy)),
460 : _mm256_mul_pd(w_c, dv_dz));
461 : __m256d rhs_v_v = _mm256_add_pd(
462 : _mm256_sub_pd(
463 : _mm256_sub_pd(
464 : _mm256_mul_pd(nu, _mm256_add_pd(_mm256_add_pd(d2v_dx2, d2v_dy2), d2v_dz2)),
465 : _mm256_mul_pd(dp_dy, rho_inv)),
466 : conv_v),
467 : src_v_v);
468 :
469 : /* RHS w: nu*(d2w_dx2+d2w_dy2+d2w_dz2) - dp_dz/rho - conv_w + src_w */
470 : __m256d conv_w = _mm256_add_pd(
471 : _mm256_add_pd(
472 : _mm256_mul_pd(u_c, dw_dx),
473 : _mm256_mul_pd(v_c, dw_dy)),
474 : _mm256_mul_pd(w_c, dw_dz));
475 : __m256d rhs_w_v = _mm256_add_pd(
476 : _mm256_sub_pd(
477 : _mm256_sub_pd(
478 : _mm256_mul_pd(nu, _mm256_add_pd(_mm256_add_pd(d2w_dx2, d2w_dy2), d2w_dz2)),
479 : _mm256_mul_pd(dp_dz, rho_inv)),
480 : conv_w),
481 : src_w_v);
482 :
483 : /* Pressure RHS: -0.1 * rho * divergence */
484 : __m256d divergence = avx2_clamp(
485 : _mm256_add_pd(_mm256_add_pd(du_dx, dv_dy), dw_dz), min_div, max_div);
486 : __m256d rhs_p_v = _mm256_mul_pd(
487 : neg_puf, _mm256_mul_pd(rho_c, divergence));
488 :
489 : /* Zero RHS for invalid lanes (rho <= eps or dx degenerate) */
490 : rhs_u_v = _mm256_and_pd(rhs_u_v, valid);
491 : rhs_v_v = _mm256_and_pd(rhs_v_v, valid);
492 : rhs_w_v = _mm256_and_pd(rhs_w_v, valid);
493 : rhs_p_v = _mm256_and_pd(rhs_p_v, valid);
494 :
495 : _mm256_storeu_pd(&rhs_u[idx], rhs_u_v);
496 : _mm256_storeu_pd(&rhs_v[idx], rhs_v_v);
497 : _mm256_storeu_pd(&rhs_w[idx], rhs_w_v);
498 : _mm256_storeu_pd(&rhs_p[idx], rhs_p_v);
499 : }
500 :
501 : /* --- Scalar: remainder from end of AVX2 to nx-2 (inclusive) ---
502 : * Handles 0-3 leftover interior points plus the periodic ir-wrap edge. */
503 : for (; i <= (ptrdiff_t)(nx - 2); i++) {
504 : size_t si = (size_t)i;
505 : size_t idx = j_off + si;
506 : size_t il = idx - 1;
507 : size_t ir = (si == nx - 2) ? j_off + 1 : idx + 1; /* periodic at last point */
508 : size_t kd = kd_off + j * nx + si;
509 : size_t ku = ku_off + j * nx + si;
510 : rk2_rhs_point(u, v, w, p, rho, T, rhs_u, rhs_v, rhs_w, rhs_p, g, params,
511 : idx, il, ir, jd_row + si, ju_row + si, kd, ku,
512 : inv_2dz, inv_dz2, si, j, k, nz, iter, dt);
513 : }
514 : }
515 : }
516 :
517 : /*
518 : * Compute RHS for all interior rows and z-planes.
519 : * k-loop is outside the OMP parallel for j-loop (MSVC OMP 2.0 compatibility).
520 : */
521 : static void compute_rhs_avx2(
522 : const double* u, const double* v, const double* w,
523 : const double* p, const double* rho, const double* T,
524 : double* rhs_u, double* rhs_v, double* rhs_w, double* rhs_p,
525 : const rk2_avx2_context_t* ctx, const grid* g,
526 : const ns_solver_params_t* params, int iter, double dt)
527 : {
528 : ptrdiff_t ny_int = (ptrdiff_t)ctx->ny;
529 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
530 : ptrdiff_t j;
531 : #ifdef _OPENMP
532 : #pragma omp parallel for schedule(static)
533 : #endif
534 : for (j = 1; j < ny_int - 1; j++) {
535 : compute_rhs_row(u, v, w, p, rho, T, rhs_u, rhs_v, rhs_w, rhs_p,
536 : ctx, g, params, (size_t)j, k, iter, dt);
537 : }
538 : }
539 : }
540 :
541 : /* ============================================================================
542 : * RK2 AVX2 MAIN LOOP
543 : * ============================================================================ */
544 :
545 : static cfd_status_t rk2_avx2_impl(flow_field* field, rk2_avx2_context_t* ctx,
546 : const grid* g, const ns_solver_params_t* params)
547 : {
548 : size_t n = ctx->nx * ctx->ny * ctx->nz;
549 : size_t bytes = n * sizeof(double);
550 : double dt = params->dt;
551 : ptrdiff_t n_int = (ptrdiff_t)n;
552 : cfd_status_t status = CFD_SUCCESS;
553 :
554 : for (int iter = 0; iter < params->max_iter; iter++) {
555 : /* Save Q^n */
556 : memcpy(ctx->u0, field->u, bytes);
557 : memcpy(ctx->v0, field->v, bytes);
558 : memcpy(ctx->w0, field->w, bytes);
559 : memcpy(ctx->p0, field->p, bytes);
560 :
561 : /* ---- Stage 1: k1 = RHS(Q^n) ---- */
562 : memset(ctx->k1_u, 0, bytes);
563 : memset(ctx->k1_v, 0, bytes);
564 : memset(ctx->k1_w, 0, bytes);
565 : memset(ctx->k1_p, 0, bytes);
566 : compute_rhs_avx2(field->u, field->v, field->w, field->p, field->rho, field->T,
567 : ctx->k1_u, ctx->k1_v, ctx->k1_w, ctx->k1_p,
568 : ctx, g, params, iter, dt);
569 :
570 : /* ---- Intermediate: field = Q^n + dt*k1 ---- */
571 : {
572 : ptrdiff_t k;
573 : #ifdef _OPENMP
574 : #pragma omp parallel for schedule(static)
575 : #endif
576 : for (k = 0; k < n_int; k++) {
577 : field->u[k] = fmax(-MAX_VELOCITY_LIMIT,
578 : fmin( MAX_VELOCITY_LIMIT,
579 : ctx->u0[k] + dt * ctx->k1_u[k]));
580 : field->v[k] = fmax(-MAX_VELOCITY_LIMIT,
581 : fmin( MAX_VELOCITY_LIMIT,
582 : ctx->v0[k] + dt * ctx->k1_v[k]));
583 : field->w[k] = fmax(-MAX_VELOCITY_LIMIT,
584 : fmin( MAX_VELOCITY_LIMIT,
585 : ctx->w0[k] + dt * ctx->k1_w[k]));
586 : field->p[k] = ctx->p0[k] + dt * ctx->k1_p[k];
587 : }
588 : }
589 :
590 : /* NOTE: Do NOT apply BCs between RK stages.
591 : * The ghost cells carry zero-derivative evolution (k1[ghost]=0),
592 : * consistent with the semi-discrete ODE. Applying BCs here would
593 : * reduce RK2 to first-order temporal accuracy. */
594 :
595 : /* ---- Stage 2: k2 = RHS(Q_pred) ---- */
596 : memset(ctx->k2_u, 0, bytes);
597 : memset(ctx->k2_v, 0, bytes);
598 : memset(ctx->k2_w, 0, bytes);
599 : memset(ctx->k2_p, 0, bytes);
600 : compute_rhs_avx2(field->u, field->v, field->w, field->p, field->rho, field->T,
601 : ctx->k2_u, ctx->k2_v, ctx->k2_w, ctx->k2_p,
602 : ctx, g, params, iter, dt);
603 :
604 : /* ---- Final update: Q^{n+1} = Q^n + (dt/2)*(k1 + k2) ---- */
605 : {
606 : double half_dt = 0.5 * dt;
607 : ptrdiff_t k;
608 : #ifdef _OPENMP
609 : #pragma omp parallel for schedule(static)
610 : #endif
611 : for (k = 0; k < n_int; k++) {
612 : field->u[k] = fmax(-MAX_VELOCITY_LIMIT,
613 : fmin( MAX_VELOCITY_LIMIT,
614 : ctx->u0[k] + half_dt * (ctx->k1_u[k] + ctx->k2_u[k])));
615 : field->v[k] = fmax(-MAX_VELOCITY_LIMIT,
616 : fmin( MAX_VELOCITY_LIMIT,
617 : ctx->v0[k] + half_dt * (ctx->k1_v[k] + ctx->k2_v[k])));
618 : field->w[k] = fmax(-MAX_VELOCITY_LIMIT,
619 : fmin( MAX_VELOCITY_LIMIT,
620 : ctx->w0[k] + half_dt * (ctx->k1_w[k] + ctx->k2_w[k])));
621 : field->p[k] = ctx->p0[k] + half_dt * (ctx->k1_p[k] + ctx->k2_p[k]);
622 : }
623 : }
624 :
625 : /* Energy equation: advance temperature after the full RK2 step */
626 : {
627 : cfd_status_t energy_status = energy_step_explicit_avx2_with_workspace(
628 : field, g, params, dt, iter * dt, ctx->T_ws, n);
629 : if (energy_status != CFD_SUCCESS) {
630 : status = energy_status;
631 : goto cleanup;
632 : }
633 : }
634 :
635 : /* Apply BCs to final state only, then configured thermal BCs */
636 : apply_boundary_conditions(field, g);
637 : status = energy_apply_thermal_bcs(field, params);
638 : if (status != CFD_SUCCESS) {
639 : goto cleanup;
640 : }
641 :
642 : /* NaN / Inf check (parallelized) */
643 : {
644 : int has_nan = 0;
645 : ptrdiff_t k;
646 : #ifdef _OPENMP
647 : #pragma omp parallel for reduction(|:has_nan) schedule(static)
648 : #endif
649 : for (k = 0; k < n_int; k++) {
650 : if (!isfinite(field->u[k]) || !isfinite(field->v[k]) ||
651 : !isfinite(field->w[k]) || !isfinite(field->p[k])) {
652 : has_nan = 1;
653 : }
654 : }
655 : if (has_nan) {
656 : status = CFD_ERROR_DIVERGED;
657 : goto cleanup;
658 : }
659 : }
660 : }
661 :
662 : cleanup:
663 : return status;
664 : }
665 :
666 : #endif /* USE_AVX2 */
667 :
668 : /* ============================================================================
669 : * PUBLIC API IMPLEMENTATIONS
670 : * ============================================================================ */
671 :
672 7 : cfd_status_t rk2_avx2_init(ns_solver_t* solver, const grid* g,
673 : const ns_solver_params_t* params)
674 : {
675 7 : (void)params;
676 :
677 : #if !USE_AVX2
678 7 : (void)solver;
679 7 : (void)g;
680 7 : cfd_set_error(CFD_ERROR_UNSUPPORTED, "AVX2 not available in this build");
681 7 : return CFD_ERROR_UNSUPPORTED;
682 : #else
683 : if (!solver || !g) {
684 : return CFD_ERROR_INVALID;
685 : }
686 : if (g->nx < 3 || g->ny < 3 || (g->nz > 1 && g->nz < 3)) {
687 : return CFD_ERROR_INVALID;
688 : }
689 :
690 : rk2_avx2_context_t* ctx =
691 : (rk2_avx2_context_t*)cfd_calloc(1, sizeof(rk2_avx2_context_t));
692 : if (!ctx) {
693 : return CFD_ERROR_NOMEM;
694 : }
695 :
696 : ctx->nx = g->nx;
697 : ctx->ny = g->ny;
698 : ctx->nz = g->nz;
699 : size_t bytes = ctx->nx * ctx->ny * g->nz * sizeof(double);
700 :
701 : /* Reject non-uniform z-spacing (solver uses constant inv_2dz/inv_dz2) */
702 : if (g->nz > 1 && g->dz) {
703 : for (size_t kk = 1; kk < g->nz - 1; kk++) {
704 : if (fabs(g->dz[kk] - g->dz[0]) > 1e-14) {
705 : cfd_free(ctx);
706 : return CFD_ERROR_INVALID;
707 : }
708 : }
709 : }
710 :
711 : /* Branch-free 3D constants */
712 : size_t plane = ctx->nx * ctx->ny;
713 : ctx->stride_z = (g->nz > 1) ? plane : 0;
714 : ctx->k_start = (g->nz > 1) ? 1 : 0;
715 : ctx->k_end = (g->nz > 1) ? (g->nz - 1) : 1;
716 : ctx->inv_2dz = (g->nz > 1 && g->dz) ? 1.0 / (2.0 * g->dz[0]) : 0.0;
717 : ctx->inv_dz2 = (g->nz > 1 && g->dz) ? 1.0 / (g->dz[0] * g->dz[0]) : 0.0;
718 :
719 : ctx->k1_u = (double*)cfd_aligned_malloc(bytes);
720 : ctx->k1_v = (double*)cfd_aligned_malloc(bytes);
721 : ctx->k1_w = (double*)cfd_aligned_malloc(bytes);
722 : ctx->k1_p = (double*)cfd_aligned_malloc(bytes);
723 : ctx->k2_u = (double*)cfd_aligned_malloc(bytes);
724 : ctx->k2_v = (double*)cfd_aligned_malloc(bytes);
725 : ctx->k2_w = (double*)cfd_aligned_malloc(bytes);
726 : ctx->k2_p = (double*)cfd_aligned_malloc(bytes);
727 : ctx->u0 = (double*)cfd_aligned_malloc(bytes);
728 : ctx->v0 = (double*)cfd_aligned_malloc(bytes);
729 : ctx->w0 = (double*)cfd_aligned_malloc(bytes);
730 : ctx->p0 = (double*)cfd_aligned_malloc(bytes);
731 : ctx->T_ws = (double*)cfd_aligned_malloc(bytes);
732 : ctx->dx_inv = (double*)cfd_aligned_malloc(ctx->nx * sizeof(double));
733 : ctx->dy_inv = (double*)cfd_aligned_malloc(ctx->ny * sizeof(double));
734 :
735 : if (!ctx->k1_u || !ctx->k1_v || !ctx->k1_w || !ctx->k1_p ||
736 : !ctx->k2_u || !ctx->k2_v || !ctx->k2_w || !ctx->k2_p ||
737 : !ctx->u0 || !ctx->v0 || !ctx->w0 || !ctx->p0 ||
738 : !ctx->T_ws || !ctx->dx_inv || !ctx->dy_inv) {
739 : cfd_aligned_free(ctx->k1_u); cfd_aligned_free(ctx->k1_v);
740 : cfd_aligned_free(ctx->k1_w); cfd_aligned_free(ctx->k1_p);
741 : cfd_aligned_free(ctx->k2_u); cfd_aligned_free(ctx->k2_v);
742 : cfd_aligned_free(ctx->k2_w); cfd_aligned_free(ctx->k2_p);
743 : cfd_aligned_free(ctx->u0); cfd_aligned_free(ctx->v0);
744 : cfd_aligned_free(ctx->w0); cfd_aligned_free(ctx->p0);
745 : cfd_aligned_free(ctx->T_ws);
746 : cfd_aligned_free(ctx->dx_inv);
747 : cfd_aligned_free(ctx->dy_inv);
748 : cfd_free(ctx);
749 : return CFD_ERROR_NOMEM;
750 : }
751 :
752 : for (size_t i = 0; i < ctx->nx; i++) {
753 : ctx->dx_inv[i] = (g->dx[i] > 1e-10) ? 1.0 / (2.0 * g->dx[i]) : 0.0;
754 : }
755 : for (size_t j = 0; j < ctx->ny; j++) {
756 : ctx->dy_inv[j] = (g->dy[j] > 1e-10) ? 1.0 / (2.0 * g->dy[j]) : 0.0;
757 : }
758 :
759 : ctx->initialized = 1;
760 : solver->context = ctx;
761 :
762 : #ifdef _OPENMP
763 : CFD_LOG_INFO("solver", "RK2 SIMD: AVX2 + OpenMP enabled (%d threads)", omp_get_max_threads());
764 : #else
765 : CFD_LOG_INFO("solver", "RK2 SIMD: AVX2 enabled (OpenMP disabled)");
766 : #endif
767 :
768 : return CFD_SUCCESS;
769 : #endif
770 : }
771 :
772 7 : void rk2_avx2_destroy(ns_solver_t* solver)
773 : {
774 7 : if (!solver || !solver->context) {
775 : return;
776 : }
777 0 : rk2_avx2_context_t* ctx = (rk2_avx2_context_t*)solver->context;
778 0 : if (ctx->initialized) {
779 0 : cfd_aligned_free(ctx->k1_u); cfd_aligned_free(ctx->k1_v);
780 0 : cfd_aligned_free(ctx->k1_w); cfd_aligned_free(ctx->k1_p);
781 0 : cfd_aligned_free(ctx->k2_u); cfd_aligned_free(ctx->k2_v);
782 0 : cfd_aligned_free(ctx->k2_w); cfd_aligned_free(ctx->k2_p);
783 0 : cfd_aligned_free(ctx->u0); cfd_aligned_free(ctx->v0);
784 0 : cfd_aligned_free(ctx->w0); cfd_aligned_free(ctx->p0);
785 0 : cfd_aligned_free(ctx->T_ws);
786 0 : cfd_aligned_free(ctx->dx_inv);
787 0 : cfd_aligned_free(ctx->dy_inv);
788 : }
789 0 : cfd_free(ctx);
790 0 : solver->context = NULL;
791 : }
792 :
793 0 : cfd_status_t rk2_avx2_step(ns_solver_t* solver, flow_field* field, const grid* g,
794 : const ns_solver_params_t* params, ns_solver_stats_t* stats)
795 : {
796 : #if !USE_AVX2
797 0 : (void)solver; (void)field; (void)g; (void)params; (void)stats;
798 0 : return CFD_ERROR_UNSUPPORTED;
799 : #else
800 : if (!solver || !solver->context || !field || !g || !params) {
801 : return CFD_ERROR_INVALID;
802 : }
803 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
804 : return CFD_ERROR_INVALID;
805 : }
806 :
807 : rk2_avx2_context_t* ctx = (rk2_avx2_context_t*)solver->context;
808 : if (field->nx != ctx->nx || field->ny != ctx->ny || field->nz != ctx->nz) {
809 : return CFD_ERROR_INVALID;
810 : }
811 :
812 : ns_solver_params_t step_params = *params;
813 : step_params.max_iter = 1;
814 :
815 : cfd_status_t status = rk2_avx2_impl(field, ctx, g, &step_params);
816 :
817 : if (stats) {
818 : stats->iterations = 1;
819 : double max_vel = 0.0, max_p = 0.0;
820 : ptrdiff_t n_s = (ptrdiff_t)(field->nx * field->ny * field->nz);
821 : double max_t = (field->T && n_s > 0) ? field->T[0] : 0.0;
822 : ptrdiff_t ks;
823 : #if defined(_OPENMP) && (_OPENMP >= 201107)
824 : #pragma omp parallel for reduction(max: max_vel, max_p, max_t) schedule(static)
825 : #endif
826 : for (ks = 0; ks < n_s; ks++) {
827 : double vel = sqrt(field->u[ks] * field->u[ks] +
828 : field->v[ks] * field->v[ks] +
829 : field->w[ks] * field->w[ks]);
830 : if (vel > max_vel) max_vel = vel;
831 : double ap = fabs(field->p[ks]);
832 : if (ap > max_p) max_p = ap;
833 : if (field->T && field->T[ks] > max_t) max_t = field->T[ks];
834 : }
835 : stats->max_velocity = max_vel;
836 : stats->max_pressure = max_p;
837 : stats->max_temperature = max_t;
838 : }
839 :
840 : return status;
841 : #endif
842 : }
843 :
844 0 : cfd_status_t rk2_avx2_solve(ns_solver_t* solver, flow_field* field, const grid* g,
845 : const ns_solver_params_t* params, ns_solver_stats_t* stats)
846 : {
847 : #if !USE_AVX2
848 0 : (void)solver; (void)field; (void)g; (void)params; (void)stats;
849 0 : return CFD_ERROR_UNSUPPORTED;
850 : #else
851 : if (!solver || !solver->context || !field || !g || !params) {
852 : return CFD_ERROR_INVALID;
853 : }
854 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
855 : return CFD_ERROR_INVALID;
856 : }
857 :
858 : rk2_avx2_context_t* ctx = (rk2_avx2_context_t*)solver->context;
859 : if (field->nx != ctx->nx || field->ny != ctx->ny || field->nz != ctx->nz) {
860 : return CFD_ERROR_INVALID;
861 : }
862 :
863 : /* Call impl directly so its internal loop runs iter=0..max_iter-1,
864 : * giving compute_source_terms the correct iteration index. */
865 : cfd_status_t status = rk2_avx2_impl(field, ctx, g, params);
866 :
867 : if (stats) {
868 : stats->iterations = params->max_iter;
869 : double max_vel = 0.0, max_p = 0.0;
870 : ptrdiff_t n_s = (ptrdiff_t)(field->nx * field->ny * field->nz);
871 : double max_t = (field->T && n_s > 0) ? field->T[0] : 0.0;
872 : ptrdiff_t ks;
873 : #if defined(_OPENMP) && (_OPENMP >= 201107)
874 : #pragma omp parallel for reduction(max: max_vel, max_p, max_t) schedule(static)
875 : #endif
876 : for (ks = 0; ks < n_s; ks++) {
877 : double vel = sqrt(field->u[ks] * field->u[ks] +
878 : field->v[ks] * field->v[ks] +
879 : field->w[ks] * field->w[ks]);
880 : if (vel > max_vel) max_vel = vel;
881 : double ap = fabs(field->p[ks]);
882 : if (ap > max_p) max_p = ap;
883 : if (field->T && field->T[ks] > max_t) max_t = field->T[ks];
884 : }
885 : stats->max_velocity = max_vel;
886 : stats->max_pressure = max_p;
887 : stats->max_temperature = max_t;
888 : }
889 : return status;
890 : #endif
891 : }
|
