Line data Source code
1 : #include "cfd/core/cfd_status.h"
2 : #include "cfd/core/grid.h"
3 : #include "cfd/core/indexing.h"
4 : #include "cfd/core/memory.h"
5 :
6 : #include "cfd/solvers/navier_stokes_solver.h"
7 : #include "cfd/solvers/energy_solver.h"
8 : #include "../../energy/energy_solver_internal.h"
9 : #include "../boundary_copy_utils.h"
10 : #include <math.h>
11 : #include <omp.h>
12 : #include <stdio.h>
13 : #include <string.h>
14 :
15 :
16 :
17 : #ifndef M_PI
18 : #define M_PI 3.14159265358979323846
19 : #endif
20 :
21 : // Physical stability limits (same as cpu solver)
22 : #define MAX_DERIVATIVE_LIMIT 100.0
23 : #define MAX_SECOND_DERIVATIVE_LIMIT 1000.0
24 : #define MAX_VELOCITY_LIMIT 100.0
25 : #define MAX_DIVERGENCE_LIMIT 10.0
26 : #define DT_CONSERVATIVE_LIMIT 0.0001
27 : #define UPDATE_LIMIT 1.0
28 : #define PRESSURE_UPDATE_FACTOR 0.1
29 :
30 : // Internal OpenMP explicit Euler implementation
31 213 : cfd_status_t explicit_euler_omp_impl(flow_field* field, const grid* grid,
32 : const ns_solver_params_t* params) {
33 213 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
34 : return CFD_ERROR_INVALID;
35 : }
36 :
37 213 : size_t nz = field->nz;
38 213 : if (nz > 1 && grid->dz) {
39 14 : for (size_t kk = 1; kk < nz - 1; kk++) {
40 12 : if (fabs(grid->dz[kk] - grid->dz[0]) > 1e-14) {
41 : return CFD_ERROR_INVALID;
42 : }
43 : }
44 : }
45 :
46 213 : size_t nx = field->nx;
47 213 : size_t ny = field->ny;
48 213 : size_t plane = nx * ny;
49 213 : size_t total = plane * nz;
50 :
51 213 : size_t stride_z = (nz > 1) ? plane : 0;
52 213 : size_t k_start = (nz > 1) ? 1 : 0;
53 213 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
54 213 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * grid->dz[0]) : 0.0;
55 2 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (grid->dz[0] * grid->dz[0]) : 0.0;
56 :
57 : // Allocate temporary arrays
58 213 : double* u_new = (double*)cfd_calloc(total, sizeof(double));
59 213 : double* v_new = (double*)cfd_calloc(total, sizeof(double));
60 213 : double* w_new = (double*)cfd_calloc(total, sizeof(double));
61 213 : double* p_new = (double*)cfd_calloc(total, sizeof(double));
62 213 : int needs_T_ws = (params->alpha > 0.0 || params->beta != 0.0);
63 213 : double* T_energy_ws = needs_T_ws
64 1 : ? (double*)cfd_calloc(total, sizeof(double)) : NULL;
65 :
66 213 : if (!u_new || !v_new || !w_new || !p_new ||
67 213 : (needs_T_ws && !T_energy_ws)) {
68 0 : cfd_free(u_new);
69 0 : cfd_free(v_new);
70 0 : cfd_free(w_new);
71 0 : cfd_free(p_new);
72 0 : cfd_free(T_energy_ws);
73 0 : return CFD_ERROR_NOMEM;
74 : }
75 :
76 : // Initialize with current values
77 213 : memcpy(u_new, field->u, total * sizeof(double));
78 213 : memcpy(v_new, field->v, total * sizeof(double));
79 213 : memcpy(w_new, field->w, total * sizeof(double));
80 213 : memcpy(p_new, field->p, total * sizeof(double));
81 :
82 213 : double conservative_dt = fmin(params->dt, DT_CONSERVATIVE_LIMIT);
83 :
84 445 : for (int iter = 0; iter < params->max_iter; iter++) {
85 474 : for (size_t kk = k_start; kk < k_end; kk++) {
86 242 : int j;
87 242 : #pragma omp parallel for schedule(static)
88 : for (j = 1; j < (int)ny - 1; j++) {
89 : for (int i = 1; i < (int)nx - 1; i++) {
90 : size_t idx = (kk * stride_z) + IDX_2D(i, j, nx);
91 :
92 : // Derivatives
93 : double du_dx = (field->u[idx + 1] - field->u[idx - 1]) / (2.0 * grid->dx[i]);
94 : double du_dy =
95 : (field->u[idx + nx] - field->u[idx - nx]) / (2.0 * grid->dy[j]);
96 : double dv_dx = (field->v[idx + 1] - field->v[idx - 1]) / (2.0 * grid->dx[i]);
97 : double dv_dy =
98 : (field->v[idx + nx] - field->v[idx - nx]) / (2.0 * grid->dy[j]);
99 :
100 : double dp_dx = (field->p[idx + 1] - field->p[idx - 1]) / (2.0 * grid->dx[i]);
101 : double dp_dy =
102 : (field->p[idx + nx] - field->p[idx - nx]) / (2.0 * grid->dy[j]);
103 :
104 : double d2u_dx2 = (field->u[idx + 1] - 2.0 * field->u[idx] + field->u[idx - 1]) /
105 : (grid->dx[i] * grid->dx[i]);
106 : double d2u_dy2 =
107 : (field->u[idx + nx] - 2.0 * field->u[idx] + field->u[idx - nx]) /
108 : (grid->dy[j] * grid->dy[j]);
109 : double d2v_dx2 = (field->v[idx + 1] - 2.0 * field->v[idx] + field->v[idx - 1]) /
110 : (grid->dx[i] * grid->dx[i]);
111 : double d2v_dy2 =
112 : (field->v[idx + nx] - 2.0 * field->v[idx] + field->v[idx - nx]) /
113 : (grid->dy[j] * grid->dy[j]);
114 :
115 : double du_dz = (field->u[idx + stride_z] - field->u[idx - stride_z]) * inv_2dz;
116 : double dv_dz = (field->v[idx + stride_z] - field->v[idx - stride_z]) * inv_2dz;
117 : double dw_dx = (field->w[idx + 1] - field->w[idx - 1]) / (2.0 * grid->dx[i]);
118 : double dw_dy = (field->w[idx + nx] - field->w[idx - nx]) / (2.0 * grid->dy[j]);
119 : double dw_dz = (field->w[idx + stride_z] - field->w[idx - stride_z]) * inv_2dz;
120 : double dp_dz = (field->p[idx + stride_z] - field->p[idx - stride_z]) * inv_2dz;
121 :
122 : double d2u_dz2 = (field->u[idx + stride_z] - 2.0 * field->u[idx] + field->u[idx - stride_z]) * inv_dz2;
123 : double d2v_dz2 = (field->v[idx + stride_z] - 2.0 * field->v[idx] + field->v[idx - stride_z]) * inv_dz2;
124 : double d2w_dx2 = (field->w[idx + 1] - 2.0 * field->w[idx] + field->w[idx - 1]) / (grid->dx[i] * grid->dx[i]);
125 : double d2w_dy2 = (field->w[idx + nx] - 2.0 * field->w[idx] + field->w[idx - nx]) / (grid->dy[j] * grid->dy[j]);
126 : double d2w_dz2 = (field->w[idx + stride_z] - 2.0 * field->w[idx] + field->w[idx - stride_z]) * inv_dz2;
127 :
128 : if (field->rho[idx] <= 1e-10) {
129 : continue;
130 : }
131 :
132 : double nu = params->mu / fmax(field->rho[idx], 1e-10);
133 : nu = fmin(nu, 1.0);
134 :
135 : // Source terms
136 : double source_u = 0.0;
137 : double source_v = 0.0;
138 : double source_w = 0.0;
139 : double z_coord = (nz > 1 && grid->z) ? grid->z[kk] : 0.0;
140 : if (params->source_func) {
141 : params->source_func(grid->x[i], grid->y[j], z_coord,
142 : iter * conservative_dt,
143 : params->source_context,
144 : &source_u, &source_v, &source_w);
145 : } else {
146 : source_u = params->source_amplitude_u * sin(M_PI * grid->y[j]) *
147 : exp(-params->source_decay_rate * iter * conservative_dt);
148 : source_v = params->source_amplitude_v * sin(2.0 * M_PI * grid->x[i]) *
149 : exp(-params->source_decay_rate * iter * conservative_dt);
150 : }
151 :
152 : // Boussinesq buoyancy source (no-op when beta == 0)
153 : energy_compute_buoyancy(field->T[idx], params,
154 : &source_u, &source_v, &source_w);
155 :
156 : // Update u
157 : double du = conservative_dt * (-field->u[idx] * du_dx - field->v[idx] * du_dy
158 : - field->w[idx] * du_dz
159 : - dp_dx / fmax(field->rho[idx], 1e-10)
160 : + nu * (d2u_dx2 + d2u_dy2 + d2u_dz2) + source_u);
161 :
162 : // Update v
163 : double dv = conservative_dt * (-field->u[idx] * dv_dx - field->v[idx] * dv_dy
164 : - field->w[idx] * dv_dz
165 : - dp_dy / fmax(field->rho[idx], 1e-10)
166 : + nu * (d2v_dx2 + d2v_dy2 + d2v_dz2) + source_v);
167 :
168 : // Update w
169 : double dw = conservative_dt * (-field->u[idx] * dw_dx - field->v[idx] * dw_dy
170 : - field->w[idx] * dw_dz
171 : - dp_dz / fmax(field->rho[idx], 1e-10)
172 : + nu * (d2w_dx2 + d2w_dy2 + d2w_dz2) + source_w);
173 :
174 : du = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, du));
175 : dv = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dv));
176 : dw = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dw));
177 :
178 : u_new[idx] = field->u[idx] + du;
179 : v_new[idx] = field->v[idx] + dv;
180 :
181 : // Limit velocity
182 : u_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, u_new[idx]));
183 : v_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, v_new[idx]));
184 : w_new[idx] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->w[idx] + dw));
185 :
186 : // Pressure update
187 : double divergence = du_dx + dv_dy + dw_dz;
188 : divergence = fmax(-MAX_DIVERGENCE_LIMIT, fmin(MAX_DIVERGENCE_LIMIT, divergence));
189 :
190 : double dp =
191 : -PRESSURE_UPDATE_FACTOR * conservative_dt * field->rho[idx] * divergence;
192 : dp = fmax(-UPDATE_LIMIT, fmin(UPDATE_LIMIT, dp));
193 : p_new[idx] = field->p[idx] + dp;
194 : }
195 : }
196 : }
197 :
198 : // Copy back (implicit barrier at end of parallel for)
199 232 : memcpy(field->u, u_new, total * sizeof(double));
200 232 : memcpy(field->v, v_new, total * sizeof(double));
201 232 : memcpy(field->w, w_new, total * sizeof(double));
202 232 : memcpy(field->p, p_new, total * sizeof(double));
203 :
204 : // Energy equation: advance temperature using updated velocity
205 : {
206 232 : cfd_status_t energy_status = energy_step_explicit_omp_with_workspace(
207 : field, grid, params, conservative_dt,
208 : iter * conservative_dt, T_energy_ws, total);
209 232 : if (energy_status != CFD_SUCCESS) {
210 0 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
211 0 : cfd_free(p_new); cfd_free(T_energy_ws);
212 0 : return energy_status;
213 : }
214 : }
215 :
216 : // Store caller-set boundary values before apply_boundary_conditions overwrites them,
217 : // then restore them afterward. Then apply configured thermal BCs.
218 232 : copy_boundary_velocities_3d(u_new, v_new, w_new,
219 232 : field->u, field->v, field->w, nx, ny, nz);
220 232 : apply_boundary_conditions(field, grid);
221 232 : copy_boundary_velocities_3d(field->u, field->v, field->w,
222 : u_new, v_new, w_new, nx, ny, nz);
223 232 : cfd_status_t bc_status = energy_apply_thermal_bcs(field, params);
224 232 : if (bc_status != CFD_SUCCESS) {
225 0 : cfd_free(u_new); cfd_free(v_new); cfd_free(w_new);
226 0 : cfd_free(p_new); cfd_free(T_energy_ws);
227 0 : return bc_status;
228 : }
229 :
230 : // Check for NaN/Inf values
231 232 : int has_nan = 0;
232 232 : ptrdiff_t total_int = (ptrdiff_t)total;
233 232 : ptrdiff_t ii;
234 232 : #pragma omp parallel for reduction(| : has_nan) schedule(static)
235 : for (ii = 0; ii < total_int; ii++) {
236 : if (!isfinite(field->u[ii]) || !isfinite(field->v[ii]) ||
237 : !isfinite(field->w[ii]) || !isfinite(field->p[ii])) {
238 : has_nan = 1;
239 : }
240 : }
241 :
242 232 : if (has_nan) {
243 0 : cfd_free(u_new);
244 0 : cfd_free(v_new);
245 0 : cfd_free(w_new);
246 0 : cfd_free(p_new);
247 0 : cfd_free(T_energy_ws);
248 0 : cfd_set_error(CFD_ERROR_DIVERGED,
249 : "NaN/Inf detected in explicit_euler_omp step");
250 0 : return CFD_ERROR_DIVERGED;
251 : }
252 : }
253 :
254 213 : cfd_free(u_new);
255 213 : cfd_free(v_new);
256 213 : cfd_free(w_new);
257 213 : cfd_free(p_new);
258 213 : cfd_free(T_energy_ws);
259 :
260 213 : return CFD_SUCCESS;
261 : }
|
