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