Line data Source code
1 : /**
2 : * Projection Method NSSolver (Chorin's Method)
3 : *
4 : * The projection method solves the incompressible Navier-Stokes equations
5 : * by splitting the solution into two steps:
6 : *
7 : * 1. Predictor Step: Compute intermediate velocity u* ignoring pressure
8 : * u* = u^n + dt * (-u·∇u + ν∇²u + f)
9 : *
10 : * 2. Pressure Projection: Solve Poisson equation for pressure
11 : * ∇²p = (ρ/dt) * ∇·u*
12 : *
13 : * 3. Corrector Step: Project velocity to be divergence-free
14 : * u^(n+1) = u* - (dt/ρ) * ∇p
15 : *
16 : * This ensures ∇·u^(n+1) = 0 (incompressibility constraint)
17 : */
18 :
19 : #include "cfd/boundary/boundary_conditions.h"
20 : #include "cfd/core/cfd_status.h"
21 : #include "cfd/core/grid.h"
22 : #include "cfd/core/indexing.h"
23 : #include "cfd/core/memory.h"
24 : #include "cfd/solvers/energy_solver.h"
25 : #include "cfd/solvers/navier_stokes_solver.h"
26 : #include "cfd/solvers/poisson_solver.h"
27 :
28 : #include "../../energy/energy_solver_internal.h"
29 : #include "../boundary_copy_utils.h"
30 :
31 : #include <math.h>
32 : #include <stdio.h>
33 : #include <string.h>
34 :
35 : #ifndef M_PI
36 : #define M_PI 3.14159265358979323846
37 : #endif
38 :
39 : // Physical limits
40 : #define MAX_VELOCITY 100.0
41 : #define MAX_PRESSURE 1000.0
42 :
43 : /**
44 : * Projection Method Solver
45 : */
46 6085 : cfd_status_t solve_projection_method(flow_field* field, const grid* grid,
47 : const ns_solver_params_t* params) {
48 6085 : if (!field || !grid || !params) {
49 : return CFD_ERROR_INVALID;
50 : }
51 6085 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
52 : return CFD_ERROR_INVALID;
53 : }
54 :
55 6085 : size_t nx = field->nx;
56 6085 : size_t ny = field->ny;
57 6085 : size_t nz = field->nz;
58 :
59 : /* Reject non-uniform z-spacing (solver uses constant dz) */
60 6085 : if (nz > 1 && grid->dz) {
61 7 : for (size_t k = 1; k < nz - 1; k++) {
62 6 : if (fabs(grid->dz[k] - grid->dz[0]) > 1e-14) {
63 : return CFD_ERROR_INVALID;
64 : }
65 : }
66 : }
67 :
68 6085 : size_t plane = nx * ny;
69 6085 : size_t total = plane * nz;
70 6085 : size_t bytes = total * sizeof(double);
71 :
72 : /* Grid spacing (uniform grid) */
73 6085 : double dx = grid->dx[0];
74 6085 : double dy = grid->dy[0];
75 6085 : double dz = (nz > 1 && grid->dz) ? grid->dz[0] : 0.0;
76 6085 : double dt = params->dt;
77 6085 : double nu = params->mu;
78 :
79 : /* Branch-free 3D constants */
80 6085 : size_t stride_z = (nz > 1) ? plane : 0;
81 6085 : size_t k_start = (nz > 1) ? 1 : 0;
82 6085 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
83 6085 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * dz) : 0.0;
84 1 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (dz * dz) : 0.0;
85 :
86 : /* Allocate temporary arrays */
87 6085 : double* u_star = (double*)cfd_calloc(total, sizeof(double));
88 6085 : double* v_star = (double*)cfd_calloc(total, sizeof(double));
89 6085 : double* w_star = (double*)cfd_calloc(total, sizeof(double));
90 6085 : double* p_new = (double*)cfd_calloc(total, sizeof(double));
91 6085 : double* p_temp = (double*)cfd_calloc(total, sizeof(double));
92 6085 : double* rhs = (double*)cfd_calloc(total, sizeof(double));
93 6085 : int needs_T_ws = (params->alpha > 0.0 || params->beta != 0.0);
94 6085 : double* T_energy_ws = needs_T_ws
95 1881 : ? (double*)cfd_calloc(total, sizeof(double)) : NULL;
96 :
97 6085 : if (!u_star || !v_star || !w_star || !p_new || !p_temp || !rhs ||
98 6085 : (needs_T_ws && !T_energy_ws)) {
99 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
100 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
101 0 : cfd_free(T_energy_ws);
102 0 : return CFD_ERROR_NOMEM;
103 : }
104 :
105 6085 : memcpy(u_star, field->u, bytes);
106 6085 : memcpy(v_star, field->v, bytes);
107 6085 : memcpy(w_star, field->w, bytes);
108 6085 : memcpy(p_new, field->p, bytes);
109 :
110 : /* Main iteration loop */
111 12208 : for (int iter = 0; iter < params->max_iter; iter++) {
112 : /* ============================================================
113 : * STEP 1: Predictor - Compute intermediate velocity u*
114 : * u* = u^n + dt * (-u·∇u + ν∇²u + f)
115 : * ============================================================ */
116 12251 : for (size_t k = k_start; k < k_end; k++) {
117 153946 : for (size_t j = 1; j < ny - 1; j++) {
118 4652628 : for (size_t i = 1; i < nx - 1; i++) {
119 4504810 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
120 :
121 4504810 : double u = field->u[idx];
122 4504810 : double v = field->v[idx];
123 4504810 : double w = field->w[idx];
124 :
125 : /* First derivatives (central) */
126 4504810 : double du_dx = (field->u[idx + 1] - field->u[idx - 1]) / (2.0 * dx);
127 4504810 : double du_dy = (field->u[idx + nx] - field->u[idx - nx]) / (2.0 * dy);
128 4504810 : double du_dz = (field->u[idx + stride_z] - field->u[idx - stride_z]) * inv_2dz;
129 :
130 4504810 : double dv_dx = (field->v[idx + 1] - field->v[idx - 1]) / (2.0 * dx);
131 4504810 : double dv_dy = (field->v[idx + nx] - field->v[idx - nx]) / (2.0 * dy);
132 4504810 : double dv_dz = (field->v[idx + stride_z] - field->v[idx - stride_z]) * inv_2dz;
133 :
134 4504810 : double dw_dx = (field->w[idx + 1] - field->w[idx - 1]) / (2.0 * dx);
135 4504810 : double dw_dy = (field->w[idx + nx] - field->w[idx - nx]) / (2.0 * dy);
136 4504810 : double dw_dz = (field->w[idx + stride_z] - field->w[idx - stride_z]) * inv_2dz;
137 :
138 : /* Convective terms: u·∇φ */
139 4504810 : double conv_u = u * du_dx + v * du_dy + w * du_dz;
140 4504810 : double conv_v = u * dv_dx + v * dv_dy + w * dv_dz;
141 4504810 : double conv_w = u * dw_dx + v * dw_dy + w * dw_dz;
142 :
143 : /* Second derivatives (viscous) */
144 4504810 : double d2u_dx2 = (field->u[idx + 1] - 2.0 * u + field->u[idx - 1]) / (dx * dx);
145 4504810 : double d2u_dy2 = (field->u[idx + nx] - 2.0 * u + field->u[idx - nx]) / (dy * dy);
146 4504810 : double d2u_dz2 = (field->u[idx + stride_z] - 2.0 * u +
147 : field->u[idx - stride_z]) * inv_dz2;
148 :
149 4504810 : double d2v_dx2 = (field->v[idx + 1] - 2.0 * v + field->v[idx - 1]) / (dx * dx);
150 4504810 : double d2v_dy2 = (field->v[idx + nx] - 2.0 * v + field->v[idx - nx]) / (dy * dy);
151 4504810 : double d2v_dz2 = (field->v[idx + stride_z] - 2.0 * v +
152 : field->v[idx - stride_z]) * inv_dz2;
153 :
154 4504810 : double d2w_dx2 = (field->w[idx + 1] - 2.0 * w + field->w[idx - 1]) / (dx * dx);
155 4504810 : double d2w_dy2 = (field->w[idx + nx] - 2.0 * w + field->w[idx - nx]) / (dy * dy);
156 4504810 : double d2w_dz2 = (field->w[idx + stride_z] - 2.0 * w +
157 : field->w[idx - stride_z]) * inv_dz2;
158 :
159 4504810 : double visc_u = nu * (d2u_dx2 + d2u_dy2 + d2u_dz2);
160 4504810 : double visc_v = nu * (d2v_dx2 + d2v_dy2 + d2v_dz2);
161 4504810 : double visc_w = nu * (d2w_dx2 + d2w_dy2 + d2w_dz2);
162 :
163 : /* Source terms */
164 4504810 : double source_u = 0.0, source_v = 0.0, source_w = 0.0;
165 4504810 : double x_coord = grid->x[i];
166 4504810 : double y_coord = grid->y[j];
167 4504810 : double z_coord = (nz > 1 && grid->z) ? grid->z[k] : 0.0;
168 4504810 : compute_source_terms(x_coord, y_coord, z_coord, iter, dt, params,
169 : &source_u, &source_v, &source_w);
170 :
171 : /* Boussinesq buoyancy source */
172 4504810 : energy_compute_buoyancy(field->T[idx], params,
173 : &source_u, &source_v, &source_w);
174 :
175 : /* Intermediate velocity (without pressure gradient) */
176 4504810 : u_star[idx] = u + dt * (-conv_u + visc_u + source_u);
177 4504810 : v_star[idx] = v + dt * (-conv_v + visc_v + source_v);
178 4504810 : w_star[idx] = w + dt * (-conv_w + visc_w + source_w);
179 :
180 4504810 : u_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, u_star[idx]));
181 4504810 : v_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, v_star[idx]));
182 4504810 : w_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, w_star[idx]));
183 : }
184 : }
185 : }
186 :
187 : /* Copy boundary values from field to star arrays */
188 6123 : copy_boundary_velocities_3d(u_star, v_star, w_star,
189 6123 : field->u, field->v, field->w, nx, ny, nz);
190 :
191 : /* ============================================================
192 : * STEP 2: Solve Poisson equation for pressure
193 : * ∇²p = (ρ/dt) * ∇·u*
194 : * ============================================================ */
195 6123 : double rho = field->rho[0];
196 6123 : if (rho < 1e-10) {
197 2 : rho = 1.0;
198 : }
199 :
200 12251 : for (size_t k = k_start; k < k_end; k++) {
201 153946 : for (size_t j = 1; j < ny - 1; j++) {
202 4652628 : for (size_t i = 1; i < nx - 1; i++) {
203 4504810 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
204 :
205 4504810 : double du_star_dx = (u_star[idx + 1] - u_star[idx - 1]) / (2.0 * dx);
206 4504810 : double dv_star_dy = (v_star[idx + nx] - v_star[idx - nx]) / (2.0 * dy);
207 4504810 : double dw_star_dz = (w_star[idx + stride_z] -
208 4504810 : w_star[idx - stride_z]) * inv_2dz;
209 :
210 4504810 : double divergence = du_star_dx + dv_star_dy + dw_star_dz;
211 4504810 : rhs[idx] = (rho / dt) * divergence;
212 : }
213 : }
214 : }
215 :
216 : /* Solve Poisson equation using library solver */
217 6123 : int poisson_iters = poisson_solve_3d(p_new, p_temp, rhs, nx, ny, nz, dx, dy, dz,
218 : POISSON_SOLVER_CG_SCALAR);
219 :
220 6123 : if (poisson_iters < 0) {
221 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
222 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
223 0 : return CFD_ERROR_MAX_ITER;
224 : }
225 :
226 : /* ============================================================
227 : * STEP 3: Corrector - Project velocity to be divergence-free
228 : * u^(n+1) = u* - (dt/ρ) * ∇p
229 : * ============================================================ */
230 6123 : double dt_over_rho = dt / rho;
231 :
232 12251 : for (size_t k = k_start; k < k_end; k++) {
233 153946 : for (size_t j = 1; j < ny - 1; j++) {
234 4652628 : for (size_t i = 1; i < nx - 1; i++) {
235 4504810 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
236 :
237 4504810 : double dp_dx = (p_new[idx + 1] - p_new[idx - 1]) / (2.0 * dx);
238 4504810 : double dp_dy = (p_new[idx + nx] - p_new[idx - nx]) / (2.0 * dy);
239 4504810 : double dp_dz = (p_new[idx + stride_z] - p_new[idx - stride_z]) * inv_2dz;
240 :
241 4504810 : field->u[idx] = u_star[idx] - dt_over_rho * dp_dx;
242 4504810 : field->v[idx] = v_star[idx] - dt_over_rho * dp_dy;
243 4504810 : field->w[idx] = w_star[idx] - dt_over_rho * dp_dz;
244 :
245 4504810 : field->u[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->u[idx]));
246 4504810 : field->v[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->v[idx]));
247 4504810 : field->w[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->w[idx]));
248 : }
249 : }
250 : }
251 :
252 : /* Update pressure */
253 6123 : memcpy(field->p, p_new, bytes);
254 :
255 : /* Energy equation: advance temperature after velocity correction */
256 : {
257 6123 : cfd_status_t energy_status = energy_step_explicit_with_workspace(
258 : field, grid, params, dt, iter * dt, T_energy_ws, total);
259 6123 : if (energy_status != CFD_SUCCESS) {
260 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
261 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
262 0 : cfd_free(T_energy_ws);
263 0 : return energy_status;
264 : }
265 : }
266 :
267 : /* Apply configured thermal BCs to temperature field */
268 6123 : cfd_status_t bc_status = energy_apply_thermal_bcs(field, params);
269 6123 : if (bc_status != CFD_SUCCESS) {
270 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
271 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
272 0 : cfd_free(T_energy_ws);
273 0 : return bc_status;
274 : }
275 :
276 : /* Restore caller-set boundary conditions */
277 6123 : copy_boundary_velocities_3d(field->u, field->v, field->w,
278 : u_star, v_star, w_star, nx, ny, nz);
279 :
280 : /* NaN/Inf check */
281 5126805 : for (size_t n = 0; n < total; n++) {
282 5120682 : if (!isfinite(field->u[n]) || !isfinite(field->v[n]) ||
283 5120682 : !isfinite(field->w[n]) || !isfinite(field->p[n])) {
284 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
285 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
286 0 : cfd_free(T_energy_ws);
287 0 : return CFD_ERROR_DIVERGED;
288 : }
289 : }
290 : }
291 :
292 6085 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
293 6085 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
294 6085 : cfd_free(T_energy_ws);
295 :
296 6085 : return CFD_SUCCESS;
297 : }
|
