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/navier_stokes_solver.h"
25 : #include "cfd/solvers/poisson_solver.h"
26 :
27 : #include "../boundary_copy_utils.h"
28 :
29 : #include <math.h>
30 : #include <stdio.h>
31 : #include <string.h>
32 :
33 : #ifndef M_PI
34 : #define M_PI 3.14159265358979323846
35 : #endif
36 :
37 : // Physical limits
38 : #define MAX_VELOCITY 100.0
39 : #define MAX_PRESSURE 1000.0
40 :
41 : /**
42 : * Projection Method Solver
43 : */
44 4753 : cfd_status_t solve_projection_method(flow_field* field, const grid* grid,
45 : const ns_solver_params_t* params) {
46 4753 : if (!field || !grid || !params) {
47 : return CFD_ERROR_INVALID;
48 : }
49 4753 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
50 : return CFD_ERROR_INVALID;
51 : }
52 :
53 4753 : size_t nx = field->nx;
54 4753 : size_t ny = field->ny;
55 4753 : size_t nz = field->nz;
56 :
57 : /* Reject non-uniform z-spacing (solver uses constant dz) */
58 4753 : if (nz > 1 && grid->dz) {
59 7857 : for (size_t k = 1; k < nz - 1; k++) {
60 7306 : if (fabs(grid->dz[k] - grid->dz[0]) > 1e-14) {
61 : return CFD_ERROR_INVALID;
62 : }
63 : }
64 : }
65 :
66 4753 : size_t plane = nx * ny;
67 4753 : size_t total = plane * nz;
68 4753 : size_t bytes = total * sizeof(double);
69 :
70 : /* Grid spacing (uniform grid) */
71 4753 : double dx = grid->dx[0];
72 4753 : double dy = grid->dy[0];
73 4753 : double dz = (nz > 1 && grid->dz) ? grid->dz[0] : 0.0;
74 4753 : double dt = params->dt;
75 4753 : double nu = params->mu;
76 :
77 : /* Branch-free 3D constants */
78 4753 : size_t stride_z = (nz > 1) ? plane : 0;
79 4753 : size_t k_start = (nz > 1) ? 1 : 0;
80 4753 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
81 4753 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * dz) : 0.0;
82 551 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (dz * dz) : 0.0;
83 :
84 : /* Allocate temporary arrays */
85 4753 : double* u_star = (double*)cfd_calloc(total, sizeof(double));
86 4753 : double* v_star = (double*)cfd_calloc(total, sizeof(double));
87 4753 : double* w_star = (double*)cfd_calloc(total, sizeof(double));
88 4753 : double* p_new = (double*)cfd_calloc(total, sizeof(double));
89 4753 : double* p_temp = (double*)cfd_calloc(total, sizeof(double));
90 4753 : double* rhs = (double*)cfd_calloc(total, sizeof(double));
91 :
92 4753 : if (!u_star || !v_star || !w_star || !p_new || !p_temp || !rhs) {
93 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
94 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
95 0 : return CFD_ERROR_NOMEM;
96 : }
97 :
98 4753 : memcpy(u_star, field->u, bytes);
99 4753 : memcpy(v_star, field->v, bytes);
100 4753 : memcpy(w_star, field->w, bytes);
101 4753 : memcpy(p_new, field->p, bytes);
102 :
103 : /* Main iteration loop */
104 9506 : for (int iter = 0; iter < params->max_iter; iter++) {
105 : /* ============================================================
106 : * STEP 1: Predictor - Compute intermediate velocity u*
107 : * u* = u^n + dt * (-u·∇u + ν∇²u + f)
108 : * ============================================================ */
109 16261 : for (size_t k = k_start; k < k_end; k++) {
110 185550 : for (size_t j = 1; j < ny - 1; j++) {
111 3201508 : for (size_t i = 1; i < nx - 1; i++) {
112 3027466 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
113 :
114 3027466 : double u = field->u[idx];
115 3027466 : double v = field->v[idx];
116 3027466 : double w = field->w[idx];
117 :
118 : /* First derivatives (central) */
119 3027466 : double du_dx = (field->u[idx + 1] - field->u[idx - 1]) / (2.0 * dx);
120 3027466 : double du_dy = (field->u[idx + nx] - field->u[idx - nx]) / (2.0 * dy);
121 3027466 : double du_dz = (field->u[idx + stride_z] - field->u[idx - stride_z]) * inv_2dz;
122 :
123 3027466 : double dv_dx = (field->v[idx + 1] - field->v[idx - 1]) / (2.0 * dx);
124 3027466 : double dv_dy = (field->v[idx + nx] - field->v[idx - nx]) / (2.0 * dy);
125 3027466 : double dv_dz = (field->v[idx + stride_z] - field->v[idx - stride_z]) * inv_2dz;
126 :
127 3027466 : double dw_dx = (field->w[idx + 1] - field->w[idx - 1]) / (2.0 * dx);
128 3027466 : double dw_dy = (field->w[idx + nx] - field->w[idx - nx]) / (2.0 * dy);
129 3027466 : double dw_dz = (field->w[idx + stride_z] - field->w[idx - stride_z]) * inv_2dz;
130 :
131 : /* Convective terms: u·∇φ */
132 3027466 : double conv_u = u * du_dx + v * du_dy + w * du_dz;
133 3027466 : double conv_v = u * dv_dx + v * dv_dy + w * dv_dz;
134 3027466 : double conv_w = u * dw_dx + v * dw_dy + w * dw_dz;
135 :
136 : /* Second derivatives (viscous) */
137 3027466 : double d2u_dx2 = (field->u[idx + 1] - 2.0 * u + field->u[idx - 1]) / (dx * dx);
138 3027466 : double d2u_dy2 = (field->u[idx + nx] - 2.0 * u + field->u[idx - nx]) / (dy * dy);
139 3027466 : double d2u_dz2 = (field->u[idx + stride_z] - 2.0 * u +
140 : field->u[idx - stride_z]) * inv_dz2;
141 :
142 3027466 : double d2v_dx2 = (field->v[idx + 1] - 2.0 * v + field->v[idx - 1]) / (dx * dx);
143 3027466 : double d2v_dy2 = (field->v[idx + nx] - 2.0 * v + field->v[idx - nx]) / (dy * dy);
144 3027466 : double d2v_dz2 = (field->v[idx + stride_z] - 2.0 * v +
145 : field->v[idx - stride_z]) * inv_dz2;
146 :
147 3027466 : double d2w_dx2 = (field->w[idx + 1] - 2.0 * w + field->w[idx - 1]) / (dx * dx);
148 3027466 : double d2w_dy2 = (field->w[idx + nx] - 2.0 * w + field->w[idx - nx]) / (dy * dy);
149 3027466 : double d2w_dz2 = (field->w[idx + stride_z] - 2.0 * w +
150 : field->w[idx - stride_z]) * inv_dz2;
151 :
152 3027466 : double visc_u = nu * (d2u_dx2 + d2u_dy2 + d2u_dz2);
153 3027466 : double visc_v = nu * (d2v_dx2 + d2v_dy2 + d2v_dz2);
154 3027466 : double visc_w = nu * (d2w_dx2 + d2w_dy2 + d2w_dz2);
155 :
156 : /* Source terms */
157 3027466 : double source_u = 0.0, source_v = 0.0, source_w = 0.0;
158 3027466 : double x_coord = grid->x[i];
159 3027466 : double y_coord = grid->y[j];
160 3027466 : double z_coord = (nz > 1 && grid->z) ? grid->z[k] : 0.0;
161 3027466 : compute_source_terms(x_coord, y_coord, z_coord, iter, dt, params,
162 : &source_u, &source_v, &source_w);
163 :
164 : /* Intermediate velocity (without pressure gradient) */
165 3027466 : u_star[idx] = u + dt * (-conv_u + visc_u + source_u);
166 3027466 : v_star[idx] = v + dt * (-conv_v + visc_v + source_v);
167 3027466 : w_star[idx] = w + dt * (-conv_w + visc_w + source_w);
168 :
169 3027466 : u_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, u_star[idx]));
170 3027466 : v_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, v_star[idx]));
171 3027466 : w_star[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, w_star[idx]));
172 : }
173 : }
174 : }
175 :
176 : /* Copy boundary values from field to star arrays */
177 4753 : copy_boundary_velocities_3d(u_star, v_star, w_star,
178 4753 : field->u, field->v, field->w, nx, ny, nz);
179 :
180 : /* ============================================================
181 : * STEP 2: Solve Poisson equation for pressure
182 : * ∇²p = (ρ/dt) * ∇·u*
183 : * ============================================================ */
184 4753 : double rho = field->rho[0];
185 4753 : if (rho < 1e-10) {
186 2 : rho = 1.0;
187 : }
188 :
189 16261 : for (size_t k = k_start; k < k_end; k++) {
190 185550 : for (size_t j = 1; j < ny - 1; j++) {
191 3201508 : for (size_t i = 1; i < nx - 1; i++) {
192 3027466 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
193 :
194 3027466 : double du_star_dx = (u_star[idx + 1] - u_star[idx - 1]) / (2.0 * dx);
195 3027466 : double dv_star_dy = (v_star[idx + nx] - v_star[idx - nx]) / (2.0 * dy);
196 3027466 : double dw_star_dz = (w_star[idx + stride_z] -
197 3027466 : w_star[idx - stride_z]) * inv_2dz;
198 :
199 3027466 : double divergence = du_star_dx + dv_star_dy + dw_star_dz;
200 3027466 : rhs[idx] = (rho / dt) * divergence;
201 : }
202 : }
203 : }
204 :
205 : /* Solve Poisson equation using library solver */
206 4753 : int poisson_iters = poisson_solve_3d(p_new, p_temp, rhs, nx, ny, nz, dx, dy, dz,
207 : POISSON_SOLVER_CG_SCALAR);
208 :
209 4753 : if (poisson_iters < 0) {
210 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
211 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
212 0 : return CFD_ERROR_MAX_ITER;
213 : }
214 :
215 : /* ============================================================
216 : * STEP 3: Corrector - Project velocity to be divergence-free
217 : * u^(n+1) = u* - (dt/ρ) * ∇p
218 : * ============================================================ */
219 4753 : double dt_over_rho = dt / rho;
220 :
221 16261 : for (size_t k = k_start; k < k_end; k++) {
222 185550 : for (size_t j = 1; j < ny - 1; j++) {
223 3201508 : for (size_t i = 1; i < nx - 1; i++) {
224 3027466 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
225 :
226 3027466 : double dp_dx = (p_new[idx + 1] - p_new[idx - 1]) / (2.0 * dx);
227 3027466 : double dp_dy = (p_new[idx + nx] - p_new[idx - nx]) / (2.0 * dy);
228 3027466 : double dp_dz = (p_new[idx + stride_z] - p_new[idx - stride_z]) * inv_2dz;
229 :
230 3027466 : field->u[idx] = u_star[idx] - dt_over_rho * dp_dx;
231 3027466 : field->v[idx] = v_star[idx] - dt_over_rho * dp_dy;
232 3027466 : field->w[idx] = w_star[idx] - dt_over_rho * dp_dz;
233 :
234 3027466 : field->u[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->u[idx]));
235 3027466 : field->v[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->v[idx]));
236 3027466 : field->w[idx] = fmax(-MAX_VELOCITY, fmin(MAX_VELOCITY, field->w[idx]));
237 : }
238 : }
239 : }
240 :
241 : /* Update pressure */
242 4753 : memcpy(field->p, p_new, bytes);
243 :
244 : /* Restore caller-set boundary conditions */
245 4753 : copy_boundary_velocities_3d(field->u, field->v, field->w,
246 : u_star, v_star, w_star, nx, ny, nz);
247 :
248 : /* NaN/Inf check */
249 4036907 : for (size_t n = 0; n < total; n++) {
250 4032154 : if (!isfinite(field->u[n]) || !isfinite(field->v[n]) ||
251 4032154 : !isfinite(field->w[n]) || !isfinite(field->p[n])) {
252 0 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
253 0 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
254 0 : return CFD_ERROR_DIVERGED;
255 : }
256 : }
257 : }
258 :
259 4753 : cfd_free(u_star); cfd_free(v_star); cfd_free(w_star);
260 4753 : cfd_free(p_new); cfd_free(p_temp); cfd_free(rhs);
261 :
262 4753 : return CFD_SUCCESS;
263 : }
|
