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