Line data Source code
1 : #include "cfd/core/cfd_status.h"
2 : #include "cfd/core/filesystem.h"
3 : #include "cfd/core/grid.h"
4 : #include "cfd/core/indexing.h"
5 : #include "cfd/core/math_utils.h"
6 : #include "cfd/core/memory.h"
7 :
8 :
9 : #include "cfd/solvers/energy_solver.h"
10 : #include "cfd/solvers/navier_stokes_solver.h"
11 :
12 : #include "../../energy/energy_solver_internal.h"
13 : #include "../boundary_copy_utils.h"
14 :
15 : #include <math.h>
16 : #include <stdio.h>
17 : #include <string.h>
18 :
19 :
20 : #ifndef M_PI
21 : #define M_PI 3.14159265358979323846
22 : #endif
23 :
24 : // Physical stability limits for numerical computation
25 : #define MAX_DERIVATIVE_LIMIT 100.0 // Maximum allowed first derivative magnitude (1/s)
26 : #define MAX_SECOND_DERIVATIVE_LIMIT 1000.0 // Maximum allowed second derivative magnitude (1/s²)
27 : #define MAX_VELOCITY_LIMIT 100.0 // Maximum allowed velocity magnitude (m/s)
28 : #define MAX_DIVERGENCE_LIMIT 10.0 // Maximum allowed velocity divergence (1/s)
29 :
30 : // Initial condition constants
31 : #define INIT_U_BASE 1.0
32 : #define INIT_U_VAR 0.1
33 : #define INIT_V_VAR 0.05
34 : #define INIT_PRESSURE 1.0
35 : #define INIT_DENSITY 1.0
36 : #define INIT_TEMP 300.0
37 :
38 : // Perturbation constants
39 : #define PERTURB_CENTER_X 1.0
40 : #define PERTURB_CENTER_Y 0.5
41 : #define PERTURB_RADIUS 0.2
42 : #define PERTURB_WIDTH_SQ 0.02
43 : #define PERTURB_MAG 0.1
44 : #define PERTURB_GRAD_FACTOR 2.0
45 :
46 : // Time stepping and stability constants
47 : #define VELOCITY_EPSILON 1e-20
48 : #define SPEED_EPSILON 1e-10
49 : #define DT_MAX_LIMIT 0.01
50 : #define DT_MIN_LIMIT 1e-6
51 : #define DT_CONSERVATIVE_LIMIT 0.0001
52 :
53 : // Update limits
54 : #define UPDATE_LIMIT 1.0
55 : #define PRESSURE_UPDATE_FACTOR 0.1
56 :
57 : // Helper function to initialize ns_solver_params_t with default values
58 239 : ns_solver_params_t ns_solver_params_default(void) {
59 239 : ns_solver_params_t params = {.dt = DEFAULT_TIME_STEP,
60 : .cfl = DEFAULT_CFL_NUMBER,
61 : .gamma = DEFAULT_GAMMA,
62 : .mu = DEFAULT_VISCOSITY,
63 : .k = DEFAULT_THERMAL_CONDUCTIVITY,
64 : .max_iter = DEFAULT_MAX_ITERATIONS,
65 : .tolerance = DEFAULT_TOLERANCE,
66 : .source_amplitude_u = DEFAULT_SOURCE_AMPLITUDE_U,
67 : .source_amplitude_v = DEFAULT_SOURCE_AMPLITUDE_V,
68 : .source_decay_rate = DEFAULT_SOURCE_DECAY_RATE,
69 : .pressure_coupling = DEFAULT_PRESSURE_COUPLING,
70 : .alpha = 0.0,
71 : .beta = 0.0,
72 : .T_ref = 0.0,
73 : .gravity = {0.0, 0.0, 0.0},
74 : .heat_source_func = NULL,
75 : .heat_source_context = NULL,
76 : .thermal_bc = {0}};
77 239 : return params;
78 : }
79 321 : flow_field* flow_field_create(size_t nx, size_t ny, size_t nz) {
80 321 : if (nx == 0 || ny == 0 || nz == 0) {
81 3 : cfd_set_error(CFD_ERROR_INVALID, "Flow field dimensions must be positive");
82 3 : return NULL;
83 : }
84 :
85 318 : flow_field* field = (flow_field*)cfd_calloc(1, sizeof(flow_field));
86 318 : if (field == NULL) {
87 : return NULL;
88 : }
89 :
90 318 : field->nx = nx;
91 318 : field->ny = ny;
92 318 : field->nz = nz;
93 :
94 318 : size_t total = nx * ny * nz;
95 :
96 : // Allocate 32-byte aligned memory for flow variables (optimized for SIMD operations)
97 318 : field->u = (double*)cfd_aligned_calloc(total, sizeof(double));
98 318 : field->v = (double*)cfd_aligned_calloc(total, sizeof(double));
99 318 : field->w = (double*)cfd_aligned_calloc(total, sizeof(double));
100 318 : field->p = (double*)cfd_aligned_calloc(total, sizeof(double));
101 318 : field->rho = (double*)cfd_aligned_calloc(total, sizeof(double));
102 318 : field->T = (double*)cfd_aligned_calloc(total, sizeof(double));
103 :
104 318 : if (!field->u || !field->v || !field->w || !field->p || !field->rho || !field->T) {
105 0 : flow_field_destroy(field);
106 0 : return NULL;
107 : }
108 :
109 : return field;
110 : }
111 :
112 319 : void flow_field_destroy(flow_field* field) {
113 319 : if (field != NULL) {
114 318 : cfd_aligned_free(field->u);
115 318 : cfd_aligned_free(field->v);
116 318 : cfd_aligned_free(field->w);
117 318 : cfd_aligned_free(field->p);
118 318 : cfd_aligned_free(field->rho);
119 318 : cfd_aligned_free(field->T);
120 318 : cfd_free(field);
121 : }
122 319 : }
123 :
124 106 : void initialize_flow_field(flow_field* field, const grid* grid) {
125 106 : size_t nx = field->nx;
126 106 : size_t ny = field->ny;
127 106 : size_t nz = field->nz;
128 106 : size_t plane = nx * ny;
129 :
130 212 : for (size_t k = 0; k < nz; k++) {
131 106 : size_t base = k * plane;
132 1097 : for (size_t j = 0; j < ny; j++) {
133 10956 : for (size_t i = 0; i < nx; i++) {
134 9965 : size_t idx = base + IDX_2D(i, j, nx);
135 9965 : double x = grid->x[i];
136 9965 : double y = grid->y[j];
137 :
138 9965 : field->u[idx] =
139 9965 : INIT_U_BASE + (INIT_U_VAR * sin(M_PI * y));
140 9965 : field->v[idx] = INIT_V_VAR * sin(2.0 * M_PI * x);
141 9965 : field->w[idx] = 0.0;
142 9965 : field->p[idx] = INIT_PRESSURE;
143 9965 : field->rho[idx] = INIT_DENSITY;
144 9965 : field->T[idx] = INIT_TEMP;
145 :
146 9965 : double cx = PERTURB_CENTER_X, cy = PERTURB_CENTER_Y;
147 9965 : double r = sqrt(((x - cx) * (x - cx)) + ((y - cy) * (y - cy)));
148 9965 : if (r < PERTURB_RADIUS) {
149 579 : field->p[idx] += PERTURB_MAG * exp(-r * r / PERTURB_WIDTH_SQ);
150 579 : double dp_dx = -PERTURB_MAG * PERTURB_GRAD_FACTOR * (x - cx) / PERTURB_WIDTH_SQ *
151 579 : exp(-r * r / PERTURB_WIDTH_SQ);
152 579 : double dp_dy = -PERTURB_MAG * PERTURB_GRAD_FACTOR * (y - cy) / PERTURB_WIDTH_SQ *
153 579 : exp(-r * r / PERTURB_WIDTH_SQ);
154 579 : field->u[idx] += -PERTURB_MAG * dp_dx;
155 579 : field->v[idx] += -PERTURB_MAG * dp_dy;
156 : }
157 : }
158 : }
159 : }
160 106 : }
161 :
162 27 : void compute_time_step(flow_field* field, const grid* grid, ns_solver_params_t* params) {
163 27 : double max_speed = 0.0;
164 27 : double dx_min = grid->dx[0];
165 27 : double dy_min = grid->dy[0];
166 :
167 : // Find minimum grid spacing
168 810 : for (size_t i = 0; i < grid->nx - 1; i++) {
169 783 : dx_min = min_double(dx_min, grid->dx[i]);
170 : }
171 810 : for (size_t j = 0; j < grid->ny - 1; j++) {
172 783 : dy_min = min_double(dy_min, grid->dy[j]);
173 : }
174 :
175 : // Find maximum wave speed
176 27 : int has_w = (grid->nz > 1 && field->w);
177 837 : for (size_t j = 0; j < field->ny; j++) {
178 35286 : for (size_t i = 0; i < field->nx; i++) {
179 34476 : size_t idx = IDX_2D(i, j, field->nx);
180 34476 : double u_speed = fabs(field->u[idx]);
181 34476 : double v_speed = fabs(field->v[idx]);
182 34476 : double sound_speed = sqrt(params->gamma * field->p[idx] / field->rho[idx]);
183 :
184 : // Optimized velocity magnitude calculation - avoid sqrt when possible
185 34476 : double vel_mag_sq = (u_speed * u_speed) + (v_speed * v_speed);
186 34476 : if (has_w) {
187 250 : double w_speed = fabs(field->w[idx]);
188 250 : vel_mag_sq += w_speed * w_speed;
189 : }
190 34476 : double vel_mag = (vel_mag_sq > VELOCITY_EPSILON) ? sqrt(vel_mag_sq) : 0.0;
191 34476 : double local_speed = vel_mag + sound_speed;
192 34476 : max_speed = max_double(max_speed, local_speed);
193 : }
194 : }
195 :
196 : // Prevent division by zero and ensure reasonable time step
197 27 : if (max_speed < SPEED_EPSILON) {
198 1 : max_speed = 1.0; // Use default if speeds are too small
199 : }
200 :
201 : // Include z-direction in CFL if 3D
202 27 : double dmin = min_double(dx_min, dy_min);
203 27 : if (grid->nz > 1 && grid->dz) {
204 4 : double dz_min = grid->dz[0];
205 60 : for (size_t k = 0; k < grid->nz - 1; k++) {
206 56 : dz_min = min_double(dz_min, grid->dz[k]);
207 : }
208 4 : dmin = min_double(dmin, dz_min);
209 : }
210 :
211 : // Compute time step based on CFL condition with safety factor
212 27 : double dt_cfl = params->cfl * dmin / max_speed;
213 :
214 : // Thermal diffusion stability constraint: dt < dx^2 / (2 * alpha * ndim)
215 27 : double dt_thermal = dt_cfl;
216 27 : if (params->alpha > 0.0) {
217 0 : int ndim = (grid->nz > 1) ? 3 : 2;
218 0 : dt_thermal = (dmin * dmin) / (2.0 * params->alpha * ndim);
219 0 : dt_thermal *= params->cfl; // Apply safety factor
220 : }
221 :
222 27 : double dt_stable = min_double(dt_cfl, dt_thermal);
223 :
224 : // Limit time step to reasonable bounds
225 27 : double dt_max = DT_MAX_LIMIT; // Maximum allowed time step
226 27 : double dt_min = DT_MIN_LIMIT; // Minimum allowed time step
227 :
228 27 : params->dt = max_double(dt_min, min_double(dt_max, dt_stable));
229 27 : }
230 :
231 11067 : void apply_boundary_conditions(flow_field* field, const grid* grid) {
232 11067 : (void)grid;
233 11067 : size_t nx = field->nx;
234 11067 : size_t ny = field->ny;
235 11067 : size_t nz = field->nz;
236 11067 : size_t plane = nx * ny;
237 :
238 : /* Apply periodic BCs in x and y for each z-plane */
239 22218 : for (size_t k = 0; k < nz; k++) {
240 11151 : size_t base = k * plane;
241 :
242 : /* x-direction periodic */
243 898246 : for (size_t j = 0; j < ny; j++) {
244 887095 : size_t left = base + IDX_2D(0, j, nx);
245 887095 : size_t right = base + IDX_2D(nx - 1, j, nx);
246 887095 : size_t src_left = base + IDX_2D(nx - 2, j, nx);
247 887095 : size_t src_right = base + IDX_2D(1, j, nx);
248 :
249 887095 : field->u[left] = field->u[src_left];
250 887095 : field->v[left] = field->v[src_left];
251 887095 : field->w[left] = field->w[src_left];
252 887095 : field->p[left] = field->p[src_left];
253 887095 : field->rho[left] = field->rho[src_left];
254 887095 : field->T[left] = field->T[src_left];
255 :
256 887095 : field->u[right] = field->u[src_right];
257 887095 : field->v[right] = field->v[src_right];
258 887095 : field->w[right] = field->w[src_right];
259 887095 : field->p[right] = field->p[src_right];
260 887095 : field->rho[right] = field->rho[src_right];
261 887095 : field->T[right] = field->T[src_right];
262 : }
263 :
264 : /* y-direction periodic */
265 901460 : for (size_t i = 0; i < nx; i++) {
266 890309 : size_t bot = base + i;
267 890309 : size_t top = base + IDX_2D(i, ny - 1, nx);
268 890309 : size_t src_bot = base + IDX_2D(i, ny - 2, nx);
269 890309 : size_t src_top = base + IDX_2D(i, 1, nx);
270 :
271 890309 : field->u[bot] = field->u[src_bot];
272 890309 : field->v[bot] = field->v[src_bot];
273 890309 : field->w[bot] = field->w[src_bot];
274 890309 : field->p[bot] = field->p[src_bot];
275 890309 : field->rho[bot] = field->rho[src_bot];
276 890309 : field->T[bot] = field->T[src_bot];
277 :
278 890309 : field->u[top] = field->u[src_top];
279 890309 : field->v[top] = field->v[src_top];
280 890309 : field->w[top] = field->w[src_top];
281 890309 : field->p[top] = field->p[src_top];
282 890309 : field->rho[top] = field->rho[src_top];
283 890309 : field->T[top] = field->T[src_top];
284 : }
285 : }
286 :
287 : /* z-direction periodic (only when nz > 1) */
288 11067 : if (nz > 1) {
289 12 : size_t front_base = 0;
290 12 : size_t back_base = (nz - 1) * plane;
291 12 : size_t src_front = (nz - 2) * plane;
292 12 : size_t src_back = 1 * plane;
293 :
294 108 : for (size_t j = 0; j < ny; j++) {
295 864 : for (size_t i = 0; i < nx; i++) {
296 768 : size_t off = IDX_2D(i, j, nx);
297 :
298 768 : field->u[front_base + off] = field->u[src_front + off];
299 768 : field->v[front_base + off] = field->v[src_front + off];
300 768 : field->w[front_base + off] = field->w[src_front + off];
301 768 : field->p[front_base + off] = field->p[src_front + off];
302 768 : field->rho[front_base + off] = field->rho[src_front + off];
303 768 : field->T[front_base + off] = field->T[src_front + off];
304 :
305 768 : field->u[back_base + off] = field->u[src_back + off];
306 768 : field->v[back_base + off] = field->v[src_back + off];
307 768 : field->w[back_base + off] = field->w[src_back + off];
308 768 : field->p[back_base + off] = field->p[src_back + off];
309 768 : field->rho[back_base + off] = field->rho[src_back + off];
310 768 : field->T[back_base + off] = field->T[src_back + off];
311 : }
312 : }
313 : }
314 11067 : }
315 :
316 : // Helper function to compute source terms consistently across all solvers
317 137695033 : void compute_source_terms(double x, double y, double z, int iter, double dt,
318 : const ns_solver_params_t* params,
319 : double* source_u, double* source_v, double* source_w) {
320 : // Use custom source function if provided
321 137695033 : if (params->source_func) {
322 126418056 : double t = iter * dt;
323 126418056 : params->source_func(x, y, z, t, params->source_context, source_u, source_v, source_w);
324 126418056 : return;
325 : }
326 :
327 : // Default source term implementation (no default w-source)
328 11276977 : *source_u =
329 11276977 : params->source_amplitude_u * sin(M_PI * y) * exp(-params->source_decay_rate * iter * dt);
330 11276977 : *source_v = params->source_amplitude_v * sin(2.0 * M_PI * x) *
331 11276977 : exp(-params->source_decay_rate * iter * dt);
332 11276977 : *source_w = 0.0;
333 : }
334 :
335 : // Internal explicit Euler implementation
336 : // This is called by the solver registry - not part of public API
337 6208 : cfd_status_t explicit_euler_impl(flow_field* field, const grid* grid, const ns_solver_params_t* params) {
338 6208 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
339 : return CFD_ERROR_INVALID;
340 : }
341 :
342 6208 : size_t nx = field->nx;
343 6208 : size_t ny = field->ny;
344 6208 : size_t nz = field->nz;
345 :
346 : /* Reject non-uniform z-spacing (solver uses constant inv_2dz/inv_dz2) */
347 6208 : if (nz > 1 && grid->dz) {
348 21 : for (size_t k = 1; k < nz - 1; k++) {
349 18 : if (fabs(grid->dz[k] - grid->dz[0]) > 1e-14) {
350 : return CFD_ERROR_INVALID;
351 : }
352 : }
353 : }
354 :
355 6208 : size_t plane = nx * ny;
356 6208 : size_t total = plane * nz;
357 6208 : size_t bytes = total * sizeof(double);
358 :
359 : /* Branch-free 3D constants: when nz==1, stride_z=0, inv_2dz=0, inv_dz2=0
360 : * so all z-terms vanish, producing identical results to 2D. */
361 6208 : size_t stride_z = (nz > 1) ? plane : 0;
362 6208 : size_t k_start = (nz > 1) ? 1 : 0;
363 6208 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
364 6208 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * grid->dz[0]) : 0.0;
365 3 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (grid->dz[0] * grid->dz[0]) : 0.0;
366 :
367 6208 : double* u_new = (double*)cfd_calloc(total, sizeof(double));
368 6208 : double* v_new = (double*)cfd_calloc(total, sizeof(double));
369 6208 : double* w_new = (double*)cfd_calloc(total, sizeof(double));
370 6208 : double* p_new = (double*)cfd_calloc(total, sizeof(double));
371 6208 : double* rho_new = (double*)cfd_calloc(total, sizeof(double));
372 6208 : int needs_T_ws = (params->alpha > 0.0 || params->beta != 0.0);
373 6208 : double* T_energy_ws = needs_T_ws
374 2 : ? (double*)cfd_calloc(total, sizeof(double)) : NULL;
375 :
376 6208 : if (!u_new || !v_new || !w_new || !p_new || !rho_new ||
377 6208 : (needs_T_ws && !T_energy_ws)) {
378 0 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
379 0 : cfd_free(p_new); cfd_free(rho_new); cfd_free(T_energy_ws);
380 0 : return CFD_ERROR_NOMEM;
381 : }
382 :
383 6208 : memcpy(u_new, field->u, bytes);
384 6208 : memcpy(v_new, field->v, bytes);
385 6208 : memcpy(w_new, field->w, bytes);
386 6208 : memcpy(p_new, field->p, bytes);
387 6208 : memcpy(rho_new, field->rho, bytes);
388 :
389 6208 : double conservative_dt = fmin(params->dt, DT_CONSERVATIVE_LIMIT);
390 :
391 12454 : for (int iter = 0; iter < params->max_iter; iter++) {
392 12507 : for (size_t k = k_start; k < k_end; k++) {
393 447384 : for (size_t j = 1; j < ny - 1; j++) {
394 45018438 : for (size_t i = 1; i < nx - 1; i++) {
395 44577315 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
396 :
397 44577315 : if (field->rho[idx] <= 1e-10) {
398 2548 : continue;
399 : }
400 44574767 : if (fabs(grid->dx[i]) < 1e-10 || fabs(grid->dy[j]) < 1e-10) {
401 0 : continue;
402 : }
403 :
404 44574767 : double u_c = field->u[idx];
405 44574767 : double v_c = field->v[idx];
406 44574767 : double w_c = field->w[idx];
407 :
408 : /* First derivatives (central differences) */
409 44574767 : double du_dx = (field->u[idx + 1] - field->u[idx - 1]) / (2.0 * grid->dx[i]);
410 44574767 : double du_dy = (field->u[idx + nx] - field->u[idx - nx]) / (2.0 * grid->dy[j]);
411 44574767 : double du_dz = (field->u[idx + stride_z] - field->u[idx - stride_z]) * inv_2dz;
412 :
413 44574767 : double dv_dx = (field->v[idx + 1] - field->v[idx - 1]) / (2.0 * grid->dx[i]);
414 44574767 : double dv_dy = (field->v[idx + nx] - field->v[idx - nx]) / (2.0 * grid->dy[j]);
415 44574767 : double dv_dz = (field->v[idx + stride_z] - field->v[idx - stride_z]) * inv_2dz;
416 :
417 44574767 : double dw_dx = (field->w[idx + 1] - field->w[idx - 1]) / (2.0 * grid->dx[i]);
418 44574767 : double dw_dy = (field->w[idx + nx] - field->w[idx - nx]) / (2.0 * grid->dy[j]);
419 44574767 : double dw_dz = (field->w[idx + stride_z] - field->w[idx - stride_z]) * inv_2dz;
420 :
421 : /* Pressure gradients */
422 44574767 : double dp_dx = (field->p[idx + 1] - field->p[idx - 1]) / (2.0 * grid->dx[i]);
423 44574767 : double dp_dy = (field->p[idx + nx] - field->p[idx - nx]) / (2.0 * grid->dy[j]);
424 44574767 : double dp_dz = (field->p[idx + stride_z] - field->p[idx - stride_z]) * inv_2dz;
425 :
426 : /* Second derivatives (viscous terms) */
427 44574767 : double d2u_dx2 = (field->u[idx + 1] - 2.0 * u_c + field->u[idx - 1]) /
428 44574767 : (grid->dx[i] * grid->dx[i]);
429 44574767 : double d2u_dy2 = (field->u[idx + nx] - 2.0 * u_c + field->u[idx - nx]) /
430 44574767 : (grid->dy[j] * grid->dy[j]);
431 44574767 : double d2u_dz2 = (field->u[idx + stride_z] - 2.0 * u_c +
432 : field->u[idx - stride_z]) * inv_dz2;
433 :
434 44574767 : double d2v_dx2 = (field->v[idx + 1] - 2.0 * v_c + field->v[idx - 1]) /
435 : (grid->dx[i] * grid->dx[i]);
436 44574767 : double d2v_dy2 = (field->v[idx + nx] - 2.0 * v_c + field->v[idx - nx]) /
437 : (grid->dy[j] * grid->dy[j]);
438 44574767 : double d2v_dz2 = (field->v[idx + stride_z] - 2.0 * v_c +
439 : field->v[idx - stride_z]) * inv_dz2;
440 :
441 44574767 : double d2w_dx2 = (field->w[idx + 1] - 2.0 * w_c + field->w[idx - 1]) /
442 : (grid->dx[i] * grid->dx[i]);
443 44574767 : double d2w_dy2 = (field->w[idx + nx] - 2.0 * w_c + field->w[idx - nx]) /
444 : (grid->dy[j] * grid->dy[j]);
445 44574767 : double d2w_dz2 = (field->w[idx + stride_z] - 2.0 * w_c +
446 : field->w[idx - stride_z]) * inv_dz2;
447 :
448 44574767 : double nu = params->mu / fmax(field->rho[idx], 1e-10);
449 44574767 : nu = fmin(nu, 1.0);
450 :
451 : /* Clamp derivatives */
452 44574767 : du_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dx));
453 44574767 : du_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dy));
454 44574767 : du_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dz));
455 44574767 : dv_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dx));
456 44574767 : dv_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dy));
457 44574767 : dv_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dz));
458 44574767 : dw_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dx));
459 44574767 : dw_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dy));
460 44574767 : dw_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dz));
461 44574767 : dp_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dx));
462 44574767 : dp_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dy));
463 44574767 : dp_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dz));
464 44574767 : d2u_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dx2));
465 44574767 : d2u_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dy2));
466 44574767 : d2u_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dz2));
467 44574767 : d2v_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dx2));
468 44574767 : d2v_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dy2));
469 44574767 : d2v_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dz2));
470 44574767 : d2w_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dx2));
471 44574767 : d2w_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dy2));
472 44574767 : d2w_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dz2));
473 :
474 : /* Source terms */
475 44574767 : double x = grid->x[i];
476 44574767 : double y = grid->y[j];
477 44574767 : double z = (nz > 1 && grid->z) ? grid->z[k] : 0.0;
478 44574767 : double source_u, source_v, source_w;
479 44574767 : compute_source_terms(x, y, z, iter, conservative_dt, params,
480 : &source_u, &source_v, &source_w);
481 :
482 : /* Boussinesq buoyancy source */
483 44574767 : energy_compute_buoyancy(field->T[idx], params,
484 : &source_u, &source_v, &source_w);
485 :
486 : /* u-momentum */
487 44574767 : double du = conservative_dt *
488 44574767 : (-u_c * du_dx - v_c * du_dy - w_c * du_dz
489 44574767 : - dp_dx / field->rho[idx]
490 44574767 : + nu * (d2u_dx2 + d2u_dy2 + d2u_dz2)
491 44574767 : + source_u);
492 :
493 : /* v-momentum */
494 44574767 : double dv = conservative_dt *
495 44574767 : (-u_c * dv_dx - v_c * dv_dy - w_c * dv_dz
496 44574767 : - dp_dy / field->rho[idx]
497 44574767 : + nu * (d2v_dx2 + d2v_dy2 + d2v_dz2)
498 44574767 : + source_v);
499 :
500 : /* w-momentum */
501 44574767 : double dw = conservative_dt *
502 44574767 : (-u_c * dw_dx - v_c * dw_dy - w_c * dw_dz
503 44574767 : - dp_dz / field->rho[idx]
504 44574767 : + nu * (d2w_dx2 + d2w_dy2 + d2w_dz2)
505 44574767 : + source_w);
506 :
507 44574767 : du = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, du));
508 44574767 : dv = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dv));
509 44574767 : dw = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dw));
510 :
511 44574767 : u_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, u_c + du));
512 44574767 : v_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, v_c + dv));
513 44574767 : w_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, w_c + dw));
514 :
515 : /* Pressure update from divergence */
516 44574767 : double divergence = du_dx + dv_dy + dw_dz;
517 44574767 : divergence = fmax(-MAX_DIVERGENCE_LIMIT, fmin(MAX_DIVERGENCE_LIMIT, divergence));
518 44574767 : double dp = -PRESSURE_UPDATE_FACTOR * conservative_dt * field->rho[idx] * divergence;
519 44574767 : dp = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dp));
520 44574767 : p_new[idx] = field->p[idx] + dp;
521 :
522 44574767 : rho_new[idx] = field->rho[idx];
523 : }
524 : }
525 : }
526 :
527 : /* Copy new solution */
528 6246 : memcpy(field->u, u_new, bytes);
529 6246 : memcpy(field->v, v_new, bytes);
530 6246 : memcpy(field->w, w_new, bytes);
531 6246 : memcpy(field->p, p_new, bytes);
532 6246 : memcpy(field->rho, rho_new, bytes);
533 :
534 : /* Energy equation: advance temperature using updated velocity */
535 : {
536 6246 : cfd_status_t energy_status = energy_step_explicit_with_workspace(
537 : field, grid, params, conservative_dt,
538 : iter * conservative_dt, T_energy_ws, total);
539 6246 : if (energy_status != CFD_SUCCESS) {
540 0 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
541 0 : cfd_free(p_new); cfd_free(rho_new); cfd_free(T_energy_ws);
542 0 : return energy_status;
543 : }
544 : }
545 :
546 : /* Save caller velocity BCs, apply periodic BCs, restore velocity BCs.
547 : * Then apply configured thermal BCs (overwrites periodic T values). */
548 6246 : copy_boundary_velocities_3d(u_new, v_new, w_new,
549 6246 : field->u, field->v, field->w, nx, ny, nz);
550 6246 : apply_boundary_conditions(field, grid);
551 6246 : copy_boundary_velocities_3d(field->u, field->v, field->w,
552 : u_new, v_new, w_new, nx, ny, nz);
553 6246 : cfd_status_t bc_status = energy_apply_thermal_bcs(field, params);
554 6246 : if (bc_status != CFD_SUCCESS) {
555 0 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
556 0 : cfd_free(p_new); cfd_free(rho_new); cfd_free(T_energy_ws);
557 0 : return bc_status;
558 : }
559 :
560 : /* NaN/Inf check */
561 46376741 : int has_nan = 0;
562 46376741 : for (size_t n = 0; n < total; n++) {
563 46370495 : if (!isfinite(field->u[n]) || !isfinite(field->v[n]) ||
564 46370495 : !isfinite(field->w[n]) || !isfinite(field->p[n])) {
565 : has_nan = 1;
566 : break;
567 : }
568 : }
569 6246 : if (has_nan) {
570 0 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
571 0 : cfd_free(p_new); cfd_free(rho_new); cfd_free(T_energy_ws);
572 0 : cfd_set_error(CFD_ERROR_DIVERGED,
573 : "NaN/Inf detected in explicit_euler step");
574 0 : return CFD_ERROR_DIVERGED;
575 : }
576 : }
577 :
578 6208 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
579 6208 : cfd_free(p_new); cfd_free(rho_new); cfd_free(T_energy_ws);
580 :
581 6208 : return CFD_SUCCESS;
582 : }
|
