Line data Source code
1 : /**
2 : * Optimized Explicit Euler NSSolver with SIMD + OpenMP
3 : *
4 : * This implementation combines SIMD vectorization (AVX2) with OpenMP
5 : * parallelization for maximum performance on multi-core CPUs.
6 : *
7 : * - Outer loops are parallelized with OpenMP
8 : * - Inner loops use AVX2 SIMD intrinsics for vectorization
9 : *
10 : * Note: When AVX2 is not enabled at compile time (CFD_ENABLE_AVX2=OFF),
11 : * this solver uses scalar code paths. Use the base explicit_euler solver
12 : * for guaranteed scalar-only execution.
13 : */
14 :
15 : // Enable C11 features for aligned_alloc
16 : #define _POSIX_C_SOURCE 200809L
17 : #define _ISOC11_SOURCE
18 : #define _USE_MATH_DEFINES
19 :
20 : #include "cfd/core/cfd_status.h"
21 : #include "cfd/core/grid.h"
22 : #include "cfd/core/indexing.h"
23 : #include "cfd/core/logging.h"
24 : #include "cfd/core/memory.h"
25 : #include "cfd/solvers/navier_stokes_solver.h"
26 : #include "cfd/solvers/energy_solver.h"
27 : #include "../../energy/energy_solver_internal.h"
28 :
29 : #include "../boundary_copy_utils.h"
30 :
31 : #include <math.h>
32 :
33 : #ifndef M_PI
34 : #define M_PI 3.14159265358979323846
35 : #endif
36 : #include <stdio.h>
37 : #include <stdlib.h>
38 : #include <string.h>
39 :
40 : #ifdef _OPENMP
41 : #include <omp.h>
42 : #endif
43 :
44 : /* AVX2 detection
45 : * CFD_HAS_AVX2 is set by CMake when -DCFD_ENABLE_AVX2=ON.
46 : * This works consistently across all compilers (GCC, Clang, MSVC).
47 : */
48 : #if defined(CFD_HAS_AVX2)
49 : #include <immintrin.h>
50 : #define USE_AVX 1
51 : #else
52 : #define USE_AVX 0
53 : #endif
54 :
55 : // Physical stability limits
56 : #define MAX_DERIVATIVE_LIMIT 100.0
57 : #define MAX_SECOND_DERIVATIVE_LIMIT 1000.0
58 : #define MAX_VELOCITY_LIMIT 100.0
59 : #define MAX_DIVERGENCE_LIMIT 10.0
60 : #define DT_CONSERVATIVE_LIMIT 0.0001
61 : #define UPDATE_LIMIT 1.0
62 : #define PRESSURE_UPDATE_FACTOR 0.1
63 :
64 : typedef struct {
65 : double* u_new;
66 : double* v_new;
67 : double* w_new;
68 : double* p_new;
69 : double* T_ws; /* Reusable scratch for the energy step (avoids per-step alloc) */
70 : double* dx_inv;
71 : double* dy_inv;
72 : size_t nx;
73 : size_t ny;
74 : size_t nz;
75 : size_t stride_z;
76 : size_t k_start;
77 : size_t k_end;
78 : double inv_2dz;
79 : double inv_dz2;
80 : int initialized;
81 : int iter_count; /* Step counter for advancing time-dependent thermal terms */
82 : } explicit_euler_simd_context;
83 :
84 : // Public API functions
85 : cfd_status_t explicit_euler_simd_init(struct NSSolver* solver, const grid* grid,
86 : const ns_solver_params_t* params);
87 : void explicit_euler_simd_destroy(struct NSSolver* solver);
88 : cfd_status_t explicit_euler_simd_step(struct NSSolver* solver, flow_field* field, const grid* grid,
89 : const ns_solver_params_t* params, ns_solver_stats_t* stats);
90 :
91 11 : cfd_status_t explicit_euler_simd_init(struct NSSolver* solver, const grid* grid,
92 : const ns_solver_params_t* params) {
93 11 : (void)params; // Unused
94 11 : if (!solver || !grid) {
95 : return CFD_ERROR_INVALID;
96 : }
97 11 : if (grid->nx < 3 || grid->ny < 3 || (grid->nz > 1 && grid->nz < 3)) {
98 : return CFD_ERROR_INVALID;
99 : }
100 :
101 11 : explicit_euler_simd_context* ctx =
102 11 : (explicit_euler_simd_context*)cfd_calloc(1, sizeof(explicit_euler_simd_context));
103 11 : if (!ctx) {
104 : return CFD_ERROR_NOMEM;
105 : }
106 :
107 11 : ctx->nx = grid->nx;
108 11 : ctx->ny = grid->ny;
109 11 : ctx->nz = grid->nz;
110 11 : size_t field_size = ctx->nx * ctx->ny * ctx->nz * sizeof(double);
111 :
112 : /* Reject non-uniform z-spacing (solver uses constant inv_2dz/inv_dz2) */
113 11 : if (grid->nz > 1 && grid->dz) {
114 14 : for (size_t kk = 1; kk < grid->nz - 1; kk++) {
115 12 : if (fabs(grid->dz[kk] - grid->dz[0]) > 1e-14) {
116 0 : cfd_free(ctx);
117 0 : return CFD_ERROR_INVALID;
118 : }
119 : }
120 : }
121 :
122 11 : size_t plane = ctx->nx * ctx->ny;
123 11 : ctx->stride_z = (grid->nz > 1) ? plane : 0;
124 11 : ctx->k_start = (grid->nz > 1) ? 1 : 0;
125 11 : ctx->k_end = (grid->nz > 1) ? (grid->nz - 1) : 1;
126 11 : ctx->inv_2dz = (grid->nz > 1 && grid->dz) ? 1.0 / (2.0 * grid->dz[0]) : 0.0;
127 11 : ctx->inv_dz2 = (grid->nz > 1 && grid->dz) ? 1.0 / (grid->dz[0] * grid->dz[0]) : 0.0;
128 :
129 11 : ctx->u_new = (double*)cfd_aligned_malloc(field_size);
130 11 : ctx->v_new = (double*)cfd_aligned_malloc(field_size);
131 11 : ctx->w_new = (double*)cfd_aligned_malloc(field_size);
132 11 : ctx->p_new = (double*)cfd_aligned_malloc(field_size);
133 11 : ctx->T_ws = (double*)cfd_aligned_malloc(field_size);
134 11 : ctx->dx_inv = (double*)cfd_aligned_malloc(ctx->nx * sizeof(double));
135 11 : ctx->dy_inv = (double*)cfd_aligned_malloc(ctx->ny * sizeof(double));
136 :
137 11 : if (!ctx->u_new || !ctx->v_new || !ctx->w_new || !ctx->p_new || !ctx->T_ws ||
138 11 : !ctx->dx_inv || !ctx->dy_inv) {
139 0 : if (ctx->u_new) {
140 0 : cfd_aligned_free(ctx->u_new);
141 : }
142 0 : if (ctx->v_new) {
143 0 : cfd_aligned_free(ctx->v_new);
144 : }
145 0 : if (ctx->w_new) {
146 0 : cfd_aligned_free(ctx->w_new);
147 : }
148 0 : if (ctx->p_new) {
149 0 : cfd_aligned_free(ctx->p_new);
150 : }
151 0 : if (ctx->T_ws) {
152 0 : cfd_aligned_free(ctx->T_ws);
153 : }
154 0 : if (ctx->dx_inv) {
155 0 : cfd_aligned_free(ctx->dx_inv);
156 : }
157 0 : if (ctx->dy_inv) {
158 0 : cfd_aligned_free(ctx->dy_inv);
159 : }
160 0 : cfd_free(ctx);
161 0 : return CFD_ERROR_NOMEM;
162 : }
163 :
164 : // Pre-compute inverses
165 273 : for (size_t i = 0; i < ctx->nx; i++) {
166 262 : ctx->dx_inv[i] = (i < ctx->nx - 1) ? 1.0 / (2.0 * grid->dx[i]) : 0.0;
167 : }
168 271 : for (size_t j = 0; j < ctx->ny; j++) {
169 260 : ctx->dy_inv[j] = (j < ctx->ny - 1) ? 1.0 / (2.0 * grid->dy[j]) : 0.0;
170 : }
171 :
172 11 : ctx->initialized = 1;
173 11 : solver->context = ctx;
174 :
175 : #if USE_AVX
176 : #ifdef _OPENMP
177 : CFD_LOG_INFO("solver", "Explicit Euler SIMD: AVX2 + OpenMP enabled (%d threads)", omp_get_max_threads());
178 : #else
179 : CFD_LOG_INFO("solver", "Explicit Euler SIMD: AVX2 enabled (OpenMP disabled)");
180 : #endif
181 : #else
182 : #ifdef _OPENMP
183 11 : CFD_LOG_INFO("solver", "Explicit Euler OMP: Scalar + OpenMP enabled (%d threads)", omp_get_max_threads());
184 : #else
185 : CFD_LOG_INFO("solver", "Explicit Euler: Scalar fallback (no SIMD or OpenMP)");
186 : #endif
187 : #endif
188 :
189 11 : return CFD_SUCCESS;
190 : }
191 :
192 15 : void explicit_euler_simd_destroy(struct NSSolver* solver) {
193 15 : if (solver && solver->context) {
194 11 : explicit_euler_simd_context* ctx = (explicit_euler_simd_context*)solver->context;
195 11 : if (ctx->initialized) {
196 11 : cfd_aligned_free(ctx->u_new);
197 11 : cfd_aligned_free(ctx->v_new);
198 11 : cfd_aligned_free(ctx->w_new);
199 11 : cfd_aligned_free(ctx->p_new);
200 11 : cfd_aligned_free(ctx->T_ws);
201 11 : cfd_aligned_free(ctx->dx_inv);
202 11 : cfd_aligned_free(ctx->dy_inv);
203 : }
204 11 : cfd_free(ctx);
205 11 : solver->context = NULL;
206 : }
207 15 : }
208 :
209 : #if USE_AVX
210 : static inline __m256d vector_fmax(__m256d a, __m256d b) {
211 : return _mm256_max_pd(a, b);
212 : }
213 : static inline __m256d vector_fmin(__m256d a, __m256d b) {
214 : return _mm256_min_pd(a, b);
215 : }
216 :
217 : typedef struct {
218 : __m256d dt_vec;
219 : __m256d max_deriv;
220 : __m256d min_deriv;
221 : __m256d max_diverg;
222 : __m256d min_diverg;
223 : __m256d max_vel_limit;
224 : __m256d min_vel_limit;
225 : __m256d one_vec;
226 : __m256d neg_one_vec;
227 : __m256d pressure_factor;
228 : __m256d two;
229 : __m256d four;
230 : __m256d epsilon;
231 : __m256d mu_vec;
232 : __m256d zero;
233 : __m256d inv_2dz_vec;
234 : __m256d inv_dz2_vec;
235 : /* Boussinesq buoyancy: accel = -beta*(T - T_ref)*g. Zero beta -> no-op. */
236 : __m256d neg_beta_vec;
237 : __m256d t_ref_vec;
238 : __m256d gx_vec;
239 : __m256d gy_vec;
240 : __m256d gz_vec;
241 : } simd_constants;
242 :
243 : static void init_simd_constants(simd_constants* c, const ns_solver_params_t* params,
244 : double conservative_dt, double inv_2dz, double inv_dz2) {
245 : c->dt_vec = _mm256_set1_pd(conservative_dt);
246 : c->max_deriv = _mm256_set1_pd(MAX_DERIVATIVE_LIMIT);
247 : c->min_deriv = _mm256_set1_pd(-MAX_DERIVATIVE_LIMIT);
248 : c->max_diverg = _mm256_set1_pd(MAX_DIVERGENCE_LIMIT);
249 : c->min_diverg = _mm256_set1_pd(-MAX_DIVERGENCE_LIMIT);
250 : c->max_vel_limit = _mm256_set1_pd(MAX_VELOCITY_LIMIT);
251 : c->min_vel_limit = _mm256_set1_pd(-MAX_VELOCITY_LIMIT);
252 : c->one_vec = _mm256_set1_pd(UPDATE_LIMIT);
253 : c->neg_one_vec = _mm256_set1_pd(-UPDATE_LIMIT);
254 : c->pressure_factor = _mm256_set1_pd(-PRESSURE_UPDATE_FACTOR);
255 : c->two = _mm256_set1_pd(2.0);
256 : c->four = _mm256_set1_pd(4.0);
257 : c->epsilon = _mm256_set1_pd(1e-10);
258 : c->mu_vec = _mm256_set1_pd(params->mu);
259 : c->zero = _mm256_setzero_pd();
260 : c->inv_2dz_vec = _mm256_set1_pd(inv_2dz);
261 : c->inv_dz2_vec = _mm256_set1_pd(inv_dz2);
262 : c->neg_beta_vec = _mm256_set1_pd(-params->beta);
263 : c->t_ref_vec = _mm256_set1_pd(params->T_ref);
264 : c->gx_vec = _mm256_set1_pd(params->gravity[0]);
265 : c->gy_vec = _mm256_set1_pd(params->gravity[1]);
266 : c->gz_vec = _mm256_set1_pd(params->gravity[2]);
267 : }
268 :
269 : static void process_simd_row(explicit_euler_simd_context* ctx, flow_field* field, const grid* grid,
270 : size_t j, const simd_constants* sc,
271 : size_t stride_z, size_t k_offset) {
272 : double dy2 = grid->dy[j] * grid->dy[j];
273 : __m256d dy_inv_val = _mm256_set1_pd(ctx->dy_inv[j]);
274 : __m256d dy2_val = _mm256_set1_pd(dy2);
275 : __m256d dy2_recip = _mm256_div_pd(sc->one_vec, dy2_val);
276 :
277 : for (size_t i = 1; i + 3 < ctx->nx - 1; i += 4) {
278 : size_t idx = k_offset + IDX_2D(i, j, ctx->nx);
279 :
280 : __m256d u = _mm256_loadu_pd(&field->u[idx]);
281 : __m256d v = _mm256_loadu_pd(&field->v[idx]);
282 : __m256d rho = _mm256_loadu_pd(&field->rho[idx]);
283 : __m256d rho_inv = _mm256_div_pd(sc->one_vec, _mm256_max_pd(rho, sc->epsilon));
284 : __m256d dx_inv_val = _mm256_loadu_pd(&ctx->dx_inv[i]);
285 :
286 : __m256d u_xp = _mm256_loadu_pd(&field->u[idx + 1]);
287 : __m256d u_xm = _mm256_loadu_pd(&field->u[idx - 1]);
288 : __m256d u_yp = _mm256_loadu_pd(&field->u[idx + ctx->nx]);
289 : __m256d u_ym = _mm256_loadu_pd(&field->u[idx - ctx->nx]);
290 : __m256d u_zp = _mm256_loadu_pd(&field->u[idx + stride_z]);
291 : __m256d u_zm = _mm256_loadu_pd(&field->u[idx - stride_z]);
292 :
293 : __m256d du_dx = _mm256_mul_pd(_mm256_sub_pd(u_xp, u_xm), dx_inv_val);
294 : __m256d du_dy = _mm256_mul_pd(_mm256_sub_pd(u_yp, u_ym), dy_inv_val);
295 : __m256d du_dz = _mm256_mul_pd(_mm256_sub_pd(u_zp, u_zm), sc->inv_2dz_vec);
296 :
297 : du_dx = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, du_dx));
298 : du_dy = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, du_dy));
299 : du_dz = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, du_dz));
300 :
301 : __m256d v_xp = _mm256_loadu_pd(&field->v[idx + 1]);
302 : __m256d v_xm = _mm256_loadu_pd(&field->v[idx - 1]);
303 : __m256d v_yp = _mm256_loadu_pd(&field->v[idx + ctx->nx]);
304 : __m256d v_ym = _mm256_loadu_pd(&field->v[idx - ctx->nx]);
305 : __m256d v_zp = _mm256_loadu_pd(&field->v[idx + stride_z]);
306 : __m256d v_zm = _mm256_loadu_pd(&field->v[idx - stride_z]);
307 :
308 : __m256d dv_dx = _mm256_mul_pd(_mm256_sub_pd(v_xp, v_xm), dx_inv_val);
309 : __m256d dv_dy = _mm256_mul_pd(_mm256_sub_pd(v_yp, v_ym), dy_inv_val);
310 : __m256d dv_dz = _mm256_mul_pd(_mm256_sub_pd(v_zp, v_zm), sc->inv_2dz_vec);
311 :
312 : dv_dx = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dv_dx));
313 : dv_dy = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dv_dy));
314 : dv_dz = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dv_dz));
315 :
316 : __m256d p_xp = _mm256_loadu_pd(&field->p[idx + 1]);
317 : __m256d p_xm = _mm256_loadu_pd(&field->p[idx - 1]);
318 : __m256d p_yp = _mm256_loadu_pd(&field->p[idx + ctx->nx]);
319 : __m256d p_ym = _mm256_loadu_pd(&field->p[idx - ctx->nx]);
320 : __m256d p_zp = _mm256_loadu_pd(&field->p[idx + stride_z]);
321 : __m256d p_zm = _mm256_loadu_pd(&field->p[idx - stride_z]);
322 :
323 : __m256d dp_dx = _mm256_mul_pd(_mm256_sub_pd(p_xp, p_xm), dx_inv_val);
324 : __m256d dp_dy = _mm256_mul_pd(_mm256_sub_pd(p_yp, p_ym), dy_inv_val);
325 : __m256d dp_dz = _mm256_mul_pd(_mm256_sub_pd(p_zp, p_zm), sc->inv_2dz_vec);
326 :
327 : dp_dx = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dp_dx));
328 : dp_dy = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dp_dy));
329 : dp_dz = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dp_dz));
330 :
331 : __m256d w = _mm256_loadu_pd(&field->w[idx]);
332 : __m256d w_xp = _mm256_loadu_pd(&field->w[idx + 1]);
333 : __m256d w_xm = _mm256_loadu_pd(&field->w[idx - 1]);
334 : __m256d w_yp = _mm256_loadu_pd(&field->w[idx + ctx->nx]);
335 : __m256d w_ym = _mm256_loadu_pd(&field->w[idx - ctx->nx]);
336 : __m256d w_zp = _mm256_loadu_pd(&field->w[idx + stride_z]);
337 : __m256d w_zm = _mm256_loadu_pd(&field->w[idx - stride_z]);
338 :
339 : __m256d dw_dx = _mm256_mul_pd(_mm256_sub_pd(w_xp, w_xm), dx_inv_val);
340 : __m256d dw_dy = _mm256_mul_pd(_mm256_sub_pd(w_yp, w_ym), dy_inv_val);
341 : __m256d dw_dz = _mm256_mul_pd(_mm256_sub_pd(w_zp, w_zm), sc->inv_2dz_vec);
342 :
343 : dw_dx = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dw_dx));
344 : dw_dy = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dw_dy));
345 : dw_dz = vector_fmax(sc->min_deriv, vector_fmin(sc->max_deriv, dw_dz));
346 :
347 : __m256d inv_dx_sq = _mm256_mul_pd(sc->four, _mm256_mul_pd(dx_inv_val, dx_inv_val));
348 :
349 : __m256d d2u_dx2 = _mm256_mul_pd(
350 : _mm256_sub_pd(_mm256_add_pd(u_xp, u_xm), _mm256_mul_pd(sc->two, u)), inv_dx_sq);
351 : __m256d d2u_dy2 = _mm256_mul_pd(
352 : _mm256_sub_pd(_mm256_add_pd(u_yp, u_ym), _mm256_mul_pd(sc->two, u)), dy2_recip);
353 : __m256d d2u_dz2 = _mm256_mul_pd(
354 : _mm256_sub_pd(_mm256_add_pd(u_zp, u_zm), _mm256_mul_pd(sc->two, u)), sc->inv_dz2_vec);
355 :
356 : __m256d d2v_dx2 = _mm256_mul_pd(
357 : _mm256_sub_pd(_mm256_add_pd(v_xp, v_xm), _mm256_mul_pd(sc->two, v)), inv_dx_sq);
358 : __m256d d2v_dy2 = _mm256_mul_pd(
359 : _mm256_sub_pd(_mm256_add_pd(v_yp, v_ym), _mm256_mul_pd(sc->two, v)), dy2_recip);
360 : __m256d d2v_dz2 = _mm256_mul_pd(
361 : _mm256_sub_pd(_mm256_add_pd(v_zp, v_zm), _mm256_mul_pd(sc->two, v)), sc->inv_dz2_vec);
362 :
363 : __m256d d2w_dx2 = _mm256_mul_pd(
364 : _mm256_sub_pd(_mm256_add_pd(w_xp, w_xm), _mm256_mul_pd(sc->two, w)), inv_dx_sq);
365 : __m256d d2w_dy2 = _mm256_mul_pd(
366 : _mm256_sub_pd(_mm256_add_pd(w_yp, w_ym), _mm256_mul_pd(sc->two, w)), dy2_recip);
367 : __m256d d2w_dz2 = _mm256_mul_pd(
368 : _mm256_sub_pd(_mm256_add_pd(w_zp, w_zm), _mm256_mul_pd(sc->two, w)), sc->inv_dz2_vec);
369 :
370 : __m256d nu = _mm256_min_pd(sc->one_vec, _mm256_mul_pd(sc->mu_vec, rho_inv));
371 :
372 : __m256d term_pres_x = _mm256_mul_pd(dp_dx, rho_inv);
373 : __m256d term_visc_u = _mm256_mul_pd(nu, _mm256_add_pd(_mm256_add_pd(d2u_dx2, d2u_dy2), d2u_dz2));
374 : __m256d conv_u = _mm256_add_pd(
375 : _mm256_add_pd(_mm256_mul_pd(u, du_dx), _mm256_mul_pd(v, du_dy)),
376 : _mm256_mul_pd(w, du_dz));
377 : __m256d du =
378 : _mm256_mul_pd(sc->dt_vec, _mm256_add_pd(_mm256_sub_pd(term_visc_u, term_pres_x),
379 : _mm256_sub_pd(sc->zero, conv_u)));
380 :
381 : __m256d term_pres_y = _mm256_mul_pd(dp_dy, rho_inv);
382 : __m256d term_visc_v = _mm256_mul_pd(nu, _mm256_add_pd(_mm256_add_pd(d2v_dx2, d2v_dy2), d2v_dz2));
383 : __m256d conv_v = _mm256_add_pd(
384 : _mm256_add_pd(_mm256_mul_pd(u, dv_dx), _mm256_mul_pd(v, dv_dy)),
385 : _mm256_mul_pd(w, dv_dz));
386 : __m256d dv =
387 : _mm256_mul_pd(sc->dt_vec, _mm256_add_pd(_mm256_sub_pd(term_visc_v, term_pres_y),
388 : _mm256_sub_pd(sc->zero, conv_v)));
389 :
390 : __m256d term_pres_z = _mm256_mul_pd(dp_dz, rho_inv);
391 : __m256d term_visc_w = _mm256_mul_pd(nu, _mm256_add_pd(_mm256_add_pd(d2w_dx2, d2w_dy2), d2w_dz2));
392 : __m256d conv_w = _mm256_add_pd(
393 : _mm256_add_pd(_mm256_mul_pd(u, dw_dx), _mm256_mul_pd(v, dw_dy)),
394 : _mm256_mul_pd(w, dw_dz));
395 : __m256d dw =
396 : _mm256_mul_pd(sc->dt_vec, _mm256_add_pd(_mm256_sub_pd(term_visc_w, term_pres_z),
397 : _mm256_sub_pd(sc->zero, conv_w)));
398 :
399 : /* Boussinesq buoyancy: dvel += dt * (-beta*(T - T_ref)) * g.
400 : * When beta == 0 this adds exactly 0.0, matching the scalar path. */
401 : __m256d Tcell = _mm256_loadu_pd(&field->T[idx]);
402 : __m256d buoy = _mm256_mul_pd(sc->neg_beta_vec, _mm256_sub_pd(Tcell, sc->t_ref_vec));
403 : __m256d dt_buoy = _mm256_mul_pd(sc->dt_vec, buoy);
404 : du = _mm256_add_pd(du, _mm256_mul_pd(dt_buoy, sc->gx_vec));
405 : dv = _mm256_add_pd(dv, _mm256_mul_pd(dt_buoy, sc->gy_vec));
406 : dw = _mm256_add_pd(dw, _mm256_mul_pd(dt_buoy, sc->gz_vec));
407 :
408 : du = vector_fmin(sc->one_vec, vector_fmax(sc->neg_one_vec, du));
409 : dv = vector_fmin(sc->one_vec, vector_fmax(sc->neg_one_vec, dv));
410 : dw = vector_fmin(sc->one_vec, vector_fmax(sc->neg_one_vec, dw));
411 :
412 : __m256d u_next = _mm256_add_pd(u, du);
413 : __m256d v_next = _mm256_add_pd(v, dv);
414 : __m256d w_next = _mm256_add_pd(w, dw);
415 :
416 : u_next = vector_fmin(sc->max_vel_limit, vector_fmax(sc->min_vel_limit, u_next));
417 : v_next = vector_fmin(sc->max_vel_limit, vector_fmax(sc->min_vel_limit, v_next));
418 : w_next = vector_fmin(sc->max_vel_limit, vector_fmax(sc->min_vel_limit, w_next));
419 :
420 : __m256d divergence = _mm256_add_pd(_mm256_add_pd(du_dx, dv_dy), dw_dz);
421 : divergence = vector_fmin(sc->max_diverg, vector_fmax(sc->min_diverg, divergence));
422 : __m256d p = _mm256_loadu_pd(&field->p[idx]);
423 : __m256d dp = _mm256_mul_pd(
424 : sc->dt_vec, _mm256_mul_pd(sc->pressure_factor, _mm256_mul_pd(rho, divergence)));
425 : dp = vector_fmin(sc->one_vec, vector_fmax(sc->neg_one_vec, dp));
426 : __m256d p_next = _mm256_add_pd(p, dp);
427 :
428 : _mm256_storeu_pd(&ctx->u_new[idx], u_next);
429 : _mm256_storeu_pd(&ctx->v_new[idx], v_next);
430 : _mm256_storeu_pd(&ctx->w_new[idx], w_next);
431 : _mm256_storeu_pd(&ctx->p_new[idx], p_next);
432 : }
433 : }
434 : #endif
435 :
436 : #if !USE_AVX
437 15676 : static void process_scalar_row(explicit_euler_simd_context* ctx, flow_field* field,
438 : const grid* grid, const ns_solver_params_t* params, size_t j,
439 : double conservative_dt, double t, size_t stride_z, size_t k_offset) {
440 877570 : for (size_t i = 1; i < ctx->nx - 1; i++) {
441 861894 : size_t idx = k_offset + IDX_2D(i, j, ctx->nx);
442 :
443 861894 : double du_dx = (field->u[idx + 1] - field->u[idx - 1]) / (2.0 * grid->dx[i]);
444 861894 : double du_dy = (field->u[idx + ctx->nx] - field->u[idx - ctx->nx]) / (2.0 * grid->dy[j]);
445 861894 : double du_dz = (field->u[idx + stride_z] - field->u[idx - stride_z]) * ctx->inv_2dz;
446 861894 : double dv_dx = (field->v[idx + 1] - field->v[idx - 1]) / (2.0 * grid->dx[i]);
447 861894 : double dv_dy = (field->v[idx + ctx->nx] - field->v[idx - ctx->nx]) / (2.0 * grid->dy[j]);
448 861894 : double dv_dz = (field->v[idx + stride_z] - field->v[idx - stride_z]) * ctx->inv_2dz;
449 861894 : double dw_dx = (field->w[idx + 1] - field->w[idx - 1]) / (2.0 * grid->dx[i]);
450 861894 : double dw_dy = (field->w[idx + ctx->nx] - field->w[idx - ctx->nx]) / (2.0 * grid->dy[j]);
451 861894 : double dw_dz = (field->w[idx + stride_z] - field->w[idx - stride_z]) * ctx->inv_2dz;
452 :
453 861894 : double dp_dx = (field->p[idx + 1] - field->p[idx - 1]) / (2.0 * grid->dx[i]);
454 861894 : double dp_dy = (field->p[idx + ctx->nx] - field->p[idx - ctx->nx]) / (2.0 * grid->dy[j]);
455 861894 : double dp_dz = (field->p[idx + stride_z] - field->p[idx - stride_z]) * ctx->inv_2dz;
456 :
457 861894 : double d2u_dx2 = (field->u[idx + 1] - 2.0 * field->u[idx] + field->u[idx - 1]) /
458 861894 : (grid->dx[i] * grid->dx[i]);
459 861894 : double d2u_dy2 = (field->u[idx + ctx->nx] - 2.0 * field->u[idx] + field->u[idx - ctx->nx]) /
460 861894 : (grid->dy[j] * grid->dy[j]);
461 861894 : double d2u_dz2 = (field->u[idx + stride_z] - 2.0 * field->u[idx] + field->u[idx - stride_z]) *
462 861894 : ctx->inv_dz2;
463 861894 : double d2v_dx2 = (field->v[idx + 1] - 2.0 * field->v[idx] + field->v[idx - 1]) /
464 : (grid->dx[i] * grid->dx[i]);
465 861894 : double d2v_dy2 = (field->v[idx + ctx->nx] - 2.0 * field->v[idx] + field->v[idx - ctx->nx]) /
466 : (grid->dy[j] * grid->dy[j]);
467 861894 : double d2v_dz2 = (field->v[idx + stride_z] - 2.0 * field->v[idx] + field->v[idx - stride_z]) *
468 : ctx->inv_dz2;
469 861894 : double d2w_dx2 = (field->w[idx + 1] - 2.0 * field->w[idx] + field->w[idx - 1]) /
470 : (grid->dx[i] * grid->dx[i]);
471 861894 : double d2w_dy2 = (field->w[idx + ctx->nx] - 2.0 * field->w[idx] + field->w[idx - ctx->nx]) /
472 : (grid->dy[j] * grid->dy[j]);
473 861894 : double d2w_dz2 = (field->w[idx + stride_z] - 2.0 * field->w[idx] + field->w[idx - stride_z]) *
474 : ctx->inv_dz2;
475 :
476 861894 : double rho = fmax(field->rho[idx], 1e-10);
477 : // Using manual define for M_PI just in case it is missed in fallback
478 : #ifndef M_PI
479 : #define M_PI 3.14159265358979323846
480 : #endif
481 861894 : double nu = fmin(params->mu / rho, 1.0);
482 :
483 861894 : du_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dx));
484 861894 : du_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dy));
485 861894 : du_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dz));
486 861894 : dv_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dx));
487 861894 : dv_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dy));
488 861894 : dv_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dz));
489 861894 : dw_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dx));
490 861894 : dw_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dy));
491 861894 : dw_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dz));
492 861894 : dp_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dx));
493 861894 : dp_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dy));
494 861894 : dp_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dz));
495 :
496 861894 : double source_u = 0.0;
497 861894 : double source_v = 0.0;
498 861894 : double source_w = 0.0;
499 861894 : if (params->source_func) {
500 0 : params->source_func(grid->x[i], grid->y[j], 0.0, t,
501 0 : params->source_context,
502 : &source_u, &source_v, &source_w);
503 861894 : } else if (params->source_amplitude_u > 0) {
504 861462 : source_u = params->source_amplitude_u * sin(M_PI * grid->y[j]);
505 861462 : source_v = params->source_amplitude_v * sin(2.0 * M_PI * grid->x[i]);
506 : }
507 :
508 : /* Boussinesq buoyancy source (no-op when beta == 0) */
509 861894 : energy_compute_buoyancy(field->T[idx], params, &source_u, &source_v, &source_w);
510 :
511 861894 : double u_c = field->u[idx];
512 861894 : double v_c = field->v[idx];
513 861894 : double w_c = field->w[idx];
514 :
515 861894 : double du = conservative_dt * (-u_c * du_dx - v_c * du_dy - w_c * du_dz -
516 861894 : dp_dx / rho + nu * (d2u_dx2 + d2u_dy2 + d2u_dz2) + source_u);
517 861894 : double dv = conservative_dt * (-u_c * dv_dx - v_c * dv_dy - w_c * dv_dz -
518 861894 : dp_dy / rho + nu * (d2v_dx2 + d2v_dy2 + d2v_dz2) + source_v);
519 861894 : double dw = conservative_dt * (-u_c * dw_dx - v_c * dw_dy - w_c * dw_dz -
520 861894 : dp_dz / rho + nu * (d2w_dx2 + d2w_dy2 + d2w_dz2) + source_w);
521 :
522 861894 : du = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, du));
523 861894 : dv = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dv));
524 861894 : dw = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dw));
525 :
526 861894 : ctx->u_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, u_c + du));
527 861894 : ctx->v_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, v_c + dv));
528 861894 : ctx->w_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, w_c + dw));
529 :
530 861894 : double divergence = du_dx + dv_dy + dw_dz;
531 861894 : divergence = fmax(-MAX_DIVERGENCE_LIMIT, fmin(MAX_DIVERGENCE_LIMIT, divergence));
532 861894 : double dp = -PRESSURE_UPDATE_FACTOR * conservative_dt * rho * divergence;
533 861894 : dp = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dp));
534 861894 : ctx->p_new[idx] = field->p[idx] + dp;
535 : }
536 15676 : }
537 : #endif
538 :
539 325 : cfd_status_t explicit_euler_simd_step(struct NSSolver* solver, flow_field* field, const grid* grid,
540 : const ns_solver_params_t* params, ns_solver_stats_t* stats) {
541 325 : if (!solver || !solver->context || !field || !grid || !params) {
542 : return CFD_ERROR_INVALID;
543 : }
544 :
545 325 : explicit_euler_simd_context* ctx = (explicit_euler_simd_context*)solver->context;
546 :
547 325 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
548 : return CFD_ERROR_INVALID;
549 : }
550 :
551 325 : if (field->nx != ctx->nx || field->ny != ctx->ny || field->nz != ctx->nz) {
552 : return CFD_ERROR_INVALID;
553 : }
554 :
555 : // Use conservative time step to match basic solver stability
556 325 : double conservative_dt = fmin(params->dt, DT_CONSERVATIVE_LIMIT);
557 :
558 : // Copy current state to temp buffers
559 325 : size_t size = ctx->nx * ctx->ny * ctx->nz;
560 325 : memcpy(ctx->u_new, field->u, size * sizeof(double));
561 325 : memcpy(ctx->v_new, field->v, size * sizeof(double));
562 325 : memcpy(ctx->w_new, field->w, size * sizeof(double));
563 325 : memcpy(ctx->p_new, field->p, size * sizeof(double));
564 :
565 325 : size_t nx = ctx->nx;
566 325 : size_t ny = ctx->ny;
567 :
568 : #if USE_AVX
569 : simd_constants sc;
570 : init_simd_constants(&sc, params, conservative_dt, ctx->inv_2dz, ctx->inv_dz2);
571 :
572 : int ny_int = (int)(ctx->ny);
573 : int j;
574 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
575 : size_t k_offset = k * ctx->stride_z;
576 : #ifdef _OPENMP
577 : #pragma omp parallel for schedule(static)
578 : #endif
579 : for (j = 1; j < ny_int - 1; j++) {
580 : process_simd_row(ctx, field, grid, (size_t)j, &sc, ctx->stride_z, k_offset);
581 : }
582 : }
583 : #else
584 325 : int ny_int = (int)(ctx->ny);
585 325 : int j;
586 660 : for (size_t k = ctx->k_start; k < ctx->k_end; k++) {
587 335 : size_t k_offset = k * ctx->stride_z;
588 : #ifdef _OPENMP
589 335 : #pragma omp parallel for schedule(static)
590 : #endif
591 : for (j = 1; j < ny_int - 1; j++) {
592 : process_scalar_row(ctx, field, grid, params, (size_t)j, conservative_dt, 0.0,
593 : ctx->stride_z, k_offset);
594 : }
595 : }
596 : #endif
597 :
598 : // Apply boundary, check NaNs, etc.
599 325 : memcpy(field->u, ctx->u_new, size * sizeof(double));
600 325 : memcpy(field->v, ctx->v_new, size * sizeof(double));
601 325 : memcpy(field->w, ctx->w_new, size * sizeof(double));
602 325 : memcpy(field->p, ctx->p_new, size * sizeof(double));
603 :
604 : // Energy equation: advance temperature using updated velocity
605 : {
606 650 : cfd_status_t energy_status = energy_step_explicit_avx2_with_workspace(
607 : field, grid, params, conservative_dt,
608 325 : ctx->iter_count * conservative_dt, ctx->T_ws, size);
609 325 : if (energy_status != CFD_SUCCESS) {
610 : return energy_status;
611 : }
612 : }
613 :
614 : // Store caller-set boundary values before apply_boundary_conditions overwrites them,
615 : // then restore them. Then apply configured thermal BCs.
616 325 : copy_boundary_velocities_3d(ctx->u_new, ctx->v_new, ctx->w_new,
617 325 : field->u, field->v, field->w, nx, ny, ctx->nz);
618 325 : apply_boundary_conditions(field, grid);
619 325 : copy_boundary_velocities_3d(field->u, field->v, field->w,
620 325 : ctx->u_new, ctx->v_new, ctx->w_new, nx, ny, ctx->nz);
621 : {
622 325 : cfd_status_t bc_status = energy_apply_thermal_bcs(field, params);
623 325 : if (bc_status != CFD_SUCCESS) {
624 : return bc_status;
625 : }
626 : }
627 :
628 325 : if (stats) {
629 305 : stats->iterations = 1;
630 : }
631 :
632 : // NaN/Inf check
633 325 : int has_nan = 0;
634 926487 : for (size_t n = 0; n < size; n++) {
635 926162 : if (!isfinite(field->u[n]) || !isfinite(field->v[n]) ||
636 926162 : !isfinite(field->w[n]) || !isfinite(field->p[n])) {
637 : has_nan = 1;
638 : break;
639 : }
640 : }
641 325 : if (has_nan) {
642 0 : cfd_set_error(CFD_ERROR_DIVERGED,
643 : "NaN/Inf detected in explicit_euler_simd step");
644 0 : return CFD_ERROR_DIVERGED;
645 : }
646 :
647 325 : ctx->iter_count++;
648 325 : return CFD_SUCCESS;
649 : }
|
