Line data Source code
1 : #include "cfd/boundary/boundary_conditions.h"
2 : #include "cfd/core/cfd_status.h"
3 : #include "cfd/core/grid.h"
4 : #include "cfd/core/indexing.h"
5 : #include "cfd/core/memory.h"
6 : #include "cfd/solvers/navier_stokes_solver.h"
7 : #include "cfd/solvers/poisson_solver.h"
8 : #include "cfd/solvers/energy_solver.h"
9 : #include "../../energy/energy_solver_internal.h"
10 :
11 : #include "../boundary_copy_utils.h"
12 :
13 : #include <math.h>
14 : #include <omp.h>
15 : #include <stdio.h>
16 : #include <string.h>
17 :
18 :
19 : #ifndef M_PI
20 : #define M_PI 3.14159265358979323846
21 : #endif
22 :
23 : // Physical limits
24 : #define MAX_VELOCITY 100.0
25 :
26 123 : cfd_status_t solve_projection_method_omp(flow_field* field, const grid* grid,
27 : const ns_solver_params_t* params) {
28 123 : if (!field || !grid || !params) {
29 : return CFD_ERROR_INVALID;
30 : }
31 123 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
32 : return CFD_ERROR_INVALID;
33 : }
34 :
35 123 : size_t nx = field->nx;
36 123 : size_t ny = field->ny;
37 123 : size_t nz = field->nz;
38 :
39 : /* Reject non-uniform z-spacing */
40 123 : if (nz > 1 && grid->dz) {
41 7 : for (size_t kk = 1; kk < nz - 1; kk++) {
42 6 : if (fabs(grid->dz[kk] - grid->dz[0]) > 1e-14) {
43 : return CFD_ERROR_INVALID;
44 : }
45 : }
46 : }
47 :
48 123 : size_t plane = nx * ny;
49 123 : size_t total = plane * nz;
50 :
51 : /* Branch-free 3D constants */
52 123 : size_t stride_z = (nz > 1) ? plane : 0;
53 123 : size_t k_start = (nz > 1) ? 1 : 0;
54 123 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
55 :
56 123 : double dx = grid->dx[0];
57 123 : double dy = grid->dy[0];
58 123 : double dz = (nz > 1 && grid->dz) ? grid->dz[0] : 0.0;
59 123 : double dt = params->dt;
60 123 : double nu = params->mu;
61 123 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * dz) : 0.0;
62 1 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (dz * dz) : 0.0;
63 :
64 123 : double* u_star = (double*)cfd_calloc(total, sizeof(double));
65 123 : double* v_star = (double*)cfd_calloc(total, sizeof(double));
66 123 : double* w_star = (double*)cfd_calloc(total, sizeof(double));
67 123 : double* p_new = (double*)cfd_calloc(total, sizeof(double));
68 123 : double* p_temp = (double*)cfd_calloc(total, sizeof(double));
69 123 : double* rhs = (double*)cfd_calloc(total, sizeof(double));
70 123 : int needs_T_ws = (params->alpha > 0.0 || params->beta != 0.0);
71 123 : double* T_energy_ws = needs_T_ws
72 1 : ? (double*)cfd_calloc(total, sizeof(double)) : NULL;
73 :
74 123 : if (!u_star || !v_star || !w_star || !p_new || !p_temp || !rhs ||
75 123 : (needs_T_ws && !T_energy_ws)) {
76 0 : cfd_free(u_star);
77 0 : cfd_free(v_star);
78 0 : cfd_free(w_star);
79 0 : cfd_free(p_new);
80 0 : cfd_free(p_temp);
81 0 : cfd_free(rhs);
82 0 : cfd_free(T_energy_ws);
83 0 : return CFD_ERROR_NOMEM;
84 : }
85 :
86 123 : memcpy(u_star, field->u, total * sizeof(double));
87 123 : memcpy(v_star, field->v, total * sizeof(double));
88 123 : memcpy(w_star, field->w, total * sizeof(double));
89 123 : memcpy(p_new, field->p, total * sizeof(double));
90 :
91 265 : for (int iter = 0; iter < params->max_iter; iter++) {
92 : /* STEP 1: Predictor — compute u_star, v_star, w_star without pressure */
93 289 : for (size_t kk = k_start; kk < k_end; kk++) {
94 147 : int j;
95 147 : #pragma omp parallel for schedule(static)
96 : for (j = 1; j < (int)ny - 1; j++) {
97 : for (int i = 1; i < (int)nx - 1; i++) {
98 : size_t idx = kk * stride_z + IDX_2D(i, j, nx);
99 :
100 : double u = field->u[idx];
101 : double v = field->v[idx];
102 : double w = field->w[idx];
103 :
104 : double du_dx = (field->u[idx + 1] - field->u[idx - 1]) / (2.0 * dx);
105 : double du_dy = (field->u[idx + nx] - field->u[idx - nx]) / (2.0 * dy);
106 : double du_dz = (field->u[idx + stride_z] - field->u[idx - stride_z]) * inv_2dz;
107 : double dv_dx = (field->v[idx + 1] - field->v[idx - 1]) / (2.0 * dx);
108 : double dv_dy = (field->v[idx + nx] - field->v[idx - nx]) / (2.0 * dy);
109 : double dv_dz = (field->v[idx + stride_z] - field->v[idx - stride_z]) * inv_2dz;
110 : double dw_dx = (field->w[idx + 1] - field->w[idx - 1]) / (2.0 * dx);
111 : double dw_dy = (field->w[idx + nx] - field->w[idx - nx]) / (2.0 * dy);
112 : double dw_dz = (field->w[idx + stride_z] - field->w[idx - stride_z]) * inv_2dz;
113 :
114 : double conv_u = u * du_dx + v * du_dy + w * du_dz;
115 : double conv_v = u * dv_dx + v * dv_dy + w * dv_dz;
116 : double conv_w = u * dw_dx + v * dw_dy + w * dw_dz;
117 :
118 : double d2u_dx2 = (field->u[idx + 1] - 2.0 * u + field->u[idx - 1]) / (dx * dx);
119 : double d2u_dy2 = (field->u[idx + nx] - 2.0 * u + field->u[idx - nx]) / (dy * dy);
120 : double d2u_dz2 = (field->u[idx + stride_z] - 2.0 * u + field->u[idx - stride_z]) * inv_dz2;
121 : double d2v_dx2 = (field->v[idx + 1] - 2.0 * v + field->v[idx - 1]) / (dx * dx);
122 : double d2v_dy2 = (field->v[idx + nx] - 2.0 * v + field->v[idx - nx]) / (dy * dy);
123 : double d2v_dz2 = (field->v[idx + stride_z] - 2.0 * v + field->v[idx - stride_z]) * inv_dz2;
124 : double d2w_dx2 = (field->w[idx + 1] - 2.0 * w + field->w[idx - 1]) / (dx * dx);
125 : double d2w_dy2 = (field->w[idx + nx] - 2.0 * w + field->w[idx - nx]) / (dy * dy);
126 : double d2w_dz2 = (field->w[idx + stride_z] - 2.0 * w + field->w[idx - stride_z]) * inv_dz2;
127 :
128 : double visc_u = nu * (d2u_dx2 + d2u_dy2 + d2u_dz2);
129 : double visc_v = nu * (d2v_dx2 + d2v_dy2 + d2v_dz2);
130 : double visc_w = nu * (d2w_dx2 + d2w_dy2 + d2w_dz2);
131 :
132 : double source_u = 0.0;
133 : double source_v = 0.0;
134 : double source_w = 0.0;
135 : double z_coord = (nz > 1 && grid->z) ? grid->z[kk] : 0.0;
136 : if (params->source_func) {
137 : params->source_func(grid->x[i], grid->y[j], z_coord, iter * dt,
138 : params->source_context,
139 : &source_u, &source_v, &source_w);
140 : } else if (params->source_amplitude_u > 0) {
141 : source_u = params->source_amplitude_u * sin(M_PI * grid->y[j]) *
142 : exp(-params->source_decay_rate * iter * dt);
143 : source_v = params->source_amplitude_v * sin(2.0 * M_PI * grid->x[i]) *
144 : exp(-params->source_decay_rate * iter * dt);
145 : }
146 :
147 : // Boussinesq buoyancy source (no-op when beta == 0)
148 : energy_compute_buoyancy(field->T[idx], params,
149 : &source_u, &source_v, &source_w);
150 :
151 : u_star[idx] = u + dt * (-conv_u + visc_u + source_u);
152 : v_star[idx] = v + dt * (-conv_v + visc_v + source_v);
153 : w_star[idx] = w + dt * (-conv_w + visc_w + source_w);
154 :
155 : u_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, u_star[idx]));
156 : v_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, v_star[idx]));
157 : w_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, w_star[idx]));
158 : }
159 : }
160 : }
161 :
162 : /* Copy boundary values from field to star arrays */
163 142 : copy_boundary_velocities_3d(u_star, v_star, w_star,
164 142 : field->u, field->v, field->w, nx, ny, nz);
165 :
166 : /* STEP 2: Pressure Poisson equation */
167 142 : double rho = field->rho[0] < 1e-10 ? 1.0 : field->rho[0];
168 :
169 289 : for (size_t kk = k_start; kk < k_end; kk++) {
170 147 : int j;
171 147 : #pragma omp parallel for schedule(static)
172 : for (j = 1; j < (int)ny - 1; j++) {
173 : for (int i = 1; i < (int)nx - 1; i++) {
174 : size_t idx = kk * stride_z + IDX_2D(i, j, nx);
175 : double du_star_dx = (u_star[idx + 1] - u_star[idx - 1]) / (2.0 * dx);
176 : double dv_star_dy = (v_star[idx + nx] - v_star[idx - nx]) / (2.0 * dy);
177 : double dw_star_dz = (w_star[idx + stride_z] - w_star[idx - stride_z]) * inv_2dz;
178 : rhs[idx] = (rho / dt) * (du_star_dx + dv_star_dy + dw_star_dz);
179 : }
180 : }
181 : }
182 :
183 : /* Use OMP CG for parallelized Poisson solve */
184 142 : int poisson_iters = poisson_solve_3d(p_new, p_temp, rhs, nx, ny, nz,
185 : dx, dy, dz, POISSON_SOLVER_CG_OMP);
186 :
187 142 : if (poisson_iters < 0) {
188 0 : cfd_free(u_star);
189 0 : cfd_free(v_star);
190 0 : cfd_free(w_star);
191 0 : cfd_free(p_new);
192 0 : cfd_free(p_temp);
193 0 : cfd_free(rhs);
194 0 : cfd_free(T_energy_ws);
195 0 : return CFD_ERROR_MAX_ITER;
196 : }
197 :
198 : /* STEP 3: Corrector — project velocities with pressure gradient */
199 142 : double dt_over_rho = dt / rho;
200 289 : for (size_t kk = k_start; kk < k_end; kk++) {
201 147 : int j;
202 147 : #pragma omp parallel for schedule(static)
203 : for (j = 1; j < (int)ny - 1; j++) {
204 : for (int i = 1; i < (int)nx - 1; i++) {
205 : size_t idx = kk * stride_z + IDX_2D(i, j, nx);
206 : double dp_dx = (p_new[idx + 1] - p_new[idx - 1]) / (2.0 * dx);
207 : double dp_dy = (p_new[idx + nx] - p_new[idx - nx]) / (2.0 * dy);
208 : double dp_dz = (p_new[idx + stride_z] - p_new[idx - stride_z]) * inv_2dz;
209 :
210 : field->u[idx] = u_star[idx] - dt_over_rho * dp_dx;
211 : field->v[idx] = v_star[idx] - dt_over_rho * dp_dy;
212 : field->w[idx] = w_star[idx] - dt_over_rho * dp_dz;
213 :
214 : field->u[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->u[idx]));
215 : field->v[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->v[idx]));
216 : field->w[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->w[idx]));
217 : }
218 : }
219 : }
220 :
221 142 : memcpy(field->p, p_new, total * sizeof(double));
222 :
223 : /* Energy equation: advance temperature after velocity correction */
224 : {
225 142 : cfd_status_t energy_status = energy_step_explicit_omp_with_workspace(
226 : field, grid, params, dt, iter * dt, T_energy_ws, total);
227 142 : if (energy_status != CFD_SUCCESS) {
228 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
229 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
230 0 : cfd_free(T_energy_ws);
231 0 : return energy_status;
232 : }
233 : }
234 :
235 : /* Apply configured thermal BCs to temperature field */
236 142 : cfd_status_t bc_status = energy_apply_thermal_bcs(field, params);
237 142 : if (bc_status != CFD_SUCCESS) {
238 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
239 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
240 0 : cfd_free(T_energy_ws);
241 0 : return bc_status;
242 : }
243 :
244 : /* Copy boundary velocity values from star arrays (which have caller's BCs) */
245 142 : copy_boundary_velocities_3d(field->u, field->v, field->w,
246 : u_star, v_star, w_star, nx, ny, nz);
247 :
248 : /* Check for NaN/Inf values (parallelized) */
249 142 : int has_nan = 0;
250 142 : ptrdiff_t total_int = (ptrdiff_t)total;
251 142 : ptrdiff_t ii;
252 142 : #pragma omp parallel for reduction(| : has_nan) schedule(static)
253 : for (ii = 0; ii < total_int; ii++) {
254 : if (!isfinite(field->u[ii]) || !isfinite(field->v[ii]) ||
255 : !isfinite(field->w[ii]) || !isfinite(field->p[ii])) {
256 : has_nan = 1;
257 : }
258 : }
259 142 : if (has_nan) {
260 0 : cfd_free(u_star);
261 0 : cfd_free(v_star);
262 0 : cfd_free(w_star);
263 0 : cfd_free(p_new);
264 0 : cfd_free(p_temp);
265 0 : cfd_free(rhs);
266 0 : cfd_free(T_energy_ws);
267 0 : return CFD_ERROR_DIVERGED;
268 : }
269 : }
270 :
271 123 : cfd_free(u_star);
272 123 : cfd_free(v_star);
273 123 : cfd_free(w_star);
274 123 : cfd_free(p_new);
275 123 : cfd_free(p_temp);
276 123 : cfd_free(rhs);
277 123 : cfd_free(T_energy_ws);
278 123 : return CFD_SUCCESS;
279 : }
|
