Line data Source code
1 : /**
2 : * @file solver_rk2.c
3 : * @brief RK2 (Heun's method) time integration for Navier-Stokes
4 : *
5 : * Second-order Runge-Kutta time stepping:
6 : * k1 = RHS(Q^n)
7 : * Q_pred = Q^n + dt * k1
8 : * k2 = RHS(Q_pred)
9 : * Q^{n+1} = Q^n + (dt/2) * (k1 + k2)
10 : *
11 : * Uses the same spatial discretisation (central differences) and physics
12 : * as the explicit Euler solver, but achieves O(dt^2) temporal accuracy
13 : * instead of O(dt).
14 : *
15 : * Branch-free 3D: when nz==1, stride_z=0 and inv_2dz/inv_dz2=0.0 cause all
16 : * z-terms to vanish, producing bit-identical results to the 2D code path.
17 : */
18 :
19 : #include "cfd/core/cfd_status.h"
20 : #include "cfd/core/grid.h"
21 : #include "cfd/core/indexing.h"
22 : #include "cfd/core/memory.h"
23 : #include "cfd/solvers/energy_solver.h"
24 : #include "cfd/solvers/navier_stokes_solver.h"
25 :
26 : #include "../../energy/energy_solver_internal.h"
27 :
28 : #include <math.h>
29 : #include <string.h>
30 :
31 : #ifndef M_PI
32 : #define M_PI 3.14159265358979323846
33 : #endif
34 :
35 : /* Physical stability limits */
36 : #define MAX_DERIVATIVE_LIMIT 100.0
37 : #define MAX_SECOND_DERIVATIVE_LIMIT 1000.0
38 : #define MAX_VELOCITY_LIMIT 100.0
39 : #define MAX_DIVERGENCE_LIMIT 10.0
40 :
41 : #define PRESSURE_UPDATE_FACTOR 0.1
42 :
43 : /* ============================================================================
44 : * RHS EVALUATION
45 : * ============================================================================ */
46 :
47 : /**
48 : * Compute the right-hand side of the semi-discrete Navier-Stokes equations.
49 : *
50 : * For each interior point, computes:
51 : * du/dt = -u·∇u - (1/ρ)∂p/∂x + ν∇²u + source_u
52 : * dv/dt = -u·∇v - (1/ρ)∂p/∂y + ν∇²v + source_v
53 : * dw/dt = -u·∇w - (1/ρ)∂p/∂z + ν∇²w + source_w
54 : * dp/dt = -0.1 * ρ * (∂u/∂x + ∂v/∂y + ∂w/∂z)
55 : *
56 : * Uses periodic stencil indices to avoid relying on ghost cells, which is
57 : * critical for preserving RK2 temporal order (ghost cell values may be stale
58 : * during intermediate RK stages).
59 : */
60 8064 : static void compute_rhs(const double* u, const double* v, const double* w,
61 : const double* p, const double* rho, const double* T,
62 : double* rhs_u, double* rhs_v, double* rhs_w, double* rhs_p,
63 : const grid* grid, const ns_solver_params_t* params,
64 : size_t nx, size_t ny, size_t nz,
65 : size_t stride_z, size_t k_start, size_t k_end,
66 : double inv_2dz, double inv_dz2,
67 : int iter, double dt) {
68 16158 : for (size_t k = k_start; k < k_end; k++) {
69 759346 : for (size_t j = 1; j < ny - 1; j++) {
70 85990684 : for (size_t i = 1; i < nx - 1; i++) {
71 85239432 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
72 :
73 : /* Safety checks */
74 85239432 : if (rho[idx] <= 1e-10) {
75 0 : rhs_u[idx] = 0.0;
76 0 : rhs_v[idx] = 0.0;
77 0 : rhs_w[idx] = 0.0;
78 0 : rhs_p[idx] = 0.0;
79 0 : continue;
80 : }
81 85239432 : if (fabs(grid->dx[i]) < 1e-10 || fabs(grid->dy[j]) < 1e-10) {
82 0 : rhs_u[idx] = 0.0;
83 0 : rhs_v[idx] = 0.0;
84 0 : rhs_w[idx] = 0.0;
85 0 : rhs_p[idx] = 0.0;
86 0 : continue;
87 : }
88 :
89 : /* Periodic stencil indices in x and y — avoids relying on ghost cells,
90 : * which is critical for preserving RK2 temporal order. */
91 85239432 : size_t il = (i > 1) ? idx - 1 : k * stride_z + IDX_2D(nx - 2, j, nx);
92 85239432 : size_t ir = (i < nx - 2) ? idx + 1 : k * stride_z + IDX_2D(1, j, nx);
93 85239432 : size_t jd = (j > 1) ? idx - nx : k * stride_z + IDX_2D(i, ny - 2, nx);
94 85239432 : size_t ju = (j < ny - 2) ? idx + nx : k * stride_z + IDX_2D(i, 1, nx);
95 :
96 : /* Periodic stencil indices in z.
97 : * When nz==1: k=0, stride_z=0, so kd=ku=idx → z-terms vanish. */
98 170478864 : size_t kd = (k > 1) ? idx - stride_z
99 85239432 : : (nz - 2) * stride_z + IDX_2D(i, j, nx);
100 170478864 : size_t ku = (k < nz - 2) ? idx + stride_z
101 85239432 : : 1 * stride_z + IDX_2D(i, j, nx);
102 :
103 : /* First derivatives (central differences) */
104 85239432 : double du_dx = (u[ir] - u[il]) / (2.0 * grid->dx[i]);
105 85239432 : double du_dy = (u[ju] - u[jd]) / (2.0 * grid->dy[j]);
106 85239432 : double du_dz = (u[ku] - u[kd]) * inv_2dz;
107 :
108 85239432 : double dv_dx = (v[ir] - v[il]) / (2.0 * grid->dx[i]);
109 85239432 : double dv_dy = (v[ju] - v[jd]) / (2.0 * grid->dy[j]);
110 85239432 : double dv_dz = (v[ku] - v[kd]) * inv_2dz;
111 :
112 85239432 : double dw_dx = (w[ir] - w[il]) / (2.0 * grid->dx[i]);
113 85239432 : double dw_dy = (w[ju] - w[jd]) / (2.0 * grid->dy[j]);
114 85239432 : double dw_dz = (w[ku] - w[kd]) * inv_2dz;
115 :
116 : /* Pressure gradients */
117 85239432 : double dp_dx = (p[ir] - p[il]) / (2.0 * grid->dx[i]);
118 85239432 : double dp_dy = (p[ju] - p[jd]) / (2.0 * grid->dy[j]);
119 85239432 : double dp_dz = (p[ku] - p[kd]) * inv_2dz;
120 :
121 : /* Second derivatives (viscous terms) */
122 85239432 : double d2u_dx2 = (u[ir] - 2.0 * u[idx] + u[il]) / (grid->dx[i] * grid->dx[i]);
123 85239432 : double d2u_dy2 = (u[ju] - 2.0 * u[idx] + u[jd]) / (grid->dy[j] * grid->dy[j]);
124 85239432 : double d2u_dz2 = (u[ku] - 2.0 * u[idx] + u[kd]) * inv_dz2;
125 :
126 85239432 : double d2v_dx2 = (v[ir] - 2.0 * v[idx] + v[il]) / (grid->dx[i] * grid->dx[i]);
127 85239432 : double d2v_dy2 = (v[ju] - 2.0 * v[idx] + v[jd]) / (grid->dy[j] * grid->dy[j]);
128 85239432 : double d2v_dz2 = (v[ku] - 2.0 * v[idx] + v[kd]) * inv_dz2;
129 :
130 85239432 : double d2w_dx2 = (w[ir] - 2.0 * w[idx] + w[il]) / (grid->dx[i] * grid->dx[i]);
131 85239432 : double d2w_dy2 = (w[ju] - 2.0 * w[idx] + w[jd]) / (grid->dy[j] * grid->dy[j]);
132 85239432 : double d2w_dz2 = (w[ku] - 2.0 * w[idx] + w[kd]) * inv_dz2;
133 :
134 : /* Kinematic viscosity */
135 85239432 : double nu = params->mu / fmax(rho[idx], 1e-10);
136 85239432 : nu = fmin(nu, 1.0);
137 :
138 : /* Clamp first derivatives */
139 85239432 : du_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dx));
140 85239432 : du_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dy));
141 85239432 : du_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, du_dz));
142 85239432 : dv_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dx));
143 85239432 : dv_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dy));
144 85239432 : dv_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dv_dz));
145 85239432 : dw_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dx));
146 85239432 : dw_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dy));
147 85239432 : dw_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dw_dz));
148 85239432 : dp_dx = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dx));
149 85239432 : dp_dy = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dy));
150 85239432 : dp_dz = fmax(-MAX_DERIVATIVE_LIMIT, fmin(MAX_DERIVATIVE_LIMIT, dp_dz));
151 :
152 : /* Clamp second derivatives */
153 85239432 : d2u_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dx2));
154 85239432 : d2u_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dy2));
155 85239432 : d2u_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2u_dz2));
156 85239432 : d2v_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dx2));
157 85239432 : d2v_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dy2));
158 85239432 : d2v_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2v_dz2));
159 85239432 : d2w_dx2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dx2));
160 85239432 : d2w_dy2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dy2));
161 85239432 : d2w_dz2 = fmax(-MAX_SECOND_DERIVATIVE_LIMIT, fmin(MAX_SECOND_DERIVATIVE_LIMIT, d2w_dz2));
162 :
163 : /* Source terms */
164 85239432 : double source_u = 0.0, source_v = 0.0, source_w = 0.0;
165 85239432 : double z_coord = (nz > 1 && grid->z) ? grid->z[k] : 0.0;
166 85239432 : compute_source_terms(grid->x[i], grid->y[j], z_coord, iter, dt,
167 : params, &source_u, &source_v, &source_w);
168 :
169 : /* Boussinesq buoyancy source */
170 85239432 : if (T) {
171 85239432 : energy_compute_buoyancy(T[idx], params,
172 : &source_u, &source_v, &source_w);
173 : }
174 :
175 : /* RHS for u-momentum */
176 85239432 : rhs_u[idx] = -u[idx] * du_dx - v[idx] * du_dy - w[idx] * du_dz
177 85239432 : - dp_dx / rho[idx]
178 85239432 : + nu * (d2u_dx2 + d2u_dy2 + d2u_dz2)
179 85239432 : + source_u;
180 :
181 : /* RHS for v-momentum */
182 85239432 : rhs_v[idx] = -u[idx] * dv_dx - v[idx] * dv_dy - w[idx] * dv_dz
183 85239432 : - dp_dy / rho[idx]
184 85239432 : + nu * (d2v_dx2 + d2v_dy2 + d2v_dz2)
185 85239432 : + source_v;
186 :
187 : /* RHS for w-momentum */
188 85239432 : rhs_w[idx] = -u[idx] * dw_dx - v[idx] * dw_dy - w[idx] * dw_dz
189 85239432 : - dp_dz / rho[idx]
190 85239432 : + nu * (d2w_dx2 + d2w_dy2 + d2w_dz2)
191 85239432 : + source_w;
192 :
193 : /* Simplified pressure RHS (divergence-based) */
194 85239432 : double divergence = du_dx + dv_dy + dw_dz;
195 85239432 : divergence = fmax(-MAX_DIVERGENCE_LIMIT,
196 : fmin(MAX_DIVERGENCE_LIMIT, divergence));
197 85239432 : rhs_p[idx] = -PRESSURE_UPDATE_FACTOR * rho[idx] * divergence;
198 : }
199 : }
200 : }
201 8064 : }
202 :
203 : /* ============================================================================
204 : * RK2 SOLVER
205 : * ============================================================================ */
206 :
207 3994 : cfd_status_t rk2_impl(flow_field* field, const grid* grid,
208 : const ns_solver_params_t* params) {
209 3994 : if (field->nx < 3 || field->ny < 3 || (field->nz > 1 && field->nz < 3)) {
210 : return CFD_ERROR_INVALID;
211 : }
212 :
213 3994 : size_t nx = field->nx;
214 3994 : size_t ny = field->ny;
215 3994 : size_t nz = field->nz;
216 :
217 : /* Reject non-uniform z-spacing (solver uses constant inv_2dz/inv_dz2) */
218 3994 : if (nz > 1 && grid->dz) {
219 21 : for (size_t k = 1; k < nz - 1; k++) {
220 18 : if (fabs(grid->dz[k] - grid->dz[0]) > 1e-14) {
221 : return CFD_ERROR_INVALID;
222 : }
223 : }
224 : }
225 :
226 3994 : size_t plane = nx * ny;
227 3994 : size_t total = plane * nz;
228 3994 : size_t bytes = total * sizeof(double);
229 :
230 : /* Branch-free 3D constants */
231 3994 : size_t stride_z = (nz > 1) ? plane : 0;
232 3994 : size_t k_start = (nz > 1) ? 1 : 0;
233 3994 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
234 3994 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * grid->dz[0]) : 0.0;
235 3 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (grid->dz[0] * grid->dz[0]) : 0.0;
236 :
237 : /* Allocate working arrays:
238 : * k1_u/v/w/p : Stage 1 derivatives
239 : * k2_u/v/w/p : Stage 2 derivatives
240 : * u0/v0/w0/p0 : Saved state Q^n
241 : */
242 3994 : double* k1_u = (double*)cfd_calloc(total, sizeof(double));
243 3994 : double* k1_v = (double*)cfd_calloc(total, sizeof(double));
244 3994 : double* k1_w = (double*)cfd_calloc(total, sizeof(double));
245 3994 : double* k1_p = (double*)cfd_calloc(total, sizeof(double));
246 3994 : double* k2_u = (double*)cfd_calloc(total, sizeof(double));
247 3994 : double* k2_v = (double*)cfd_calloc(total, sizeof(double));
248 3994 : double* k2_w = (double*)cfd_calloc(total, sizeof(double));
249 3994 : double* k2_p = (double*)cfd_calloc(total, sizeof(double));
250 3994 : double* u0 = (double*)cfd_calloc(total, sizeof(double));
251 3994 : double* v0 = (double*)cfd_calloc(total, sizeof(double));
252 3994 : double* w0 = (double*)cfd_calloc(total, sizeof(double));
253 3994 : double* p0 = (double*)cfd_calloc(total, sizeof(double));
254 3994 : int needs_T_ws = (params->alpha > 0.0 || params->beta != 0.0);
255 3994 : double* T_energy_ws = needs_T_ws
256 2 : ? (double*)cfd_calloc(total, sizeof(double)) : NULL;
257 :
258 3994 : if (!k1_u || !k1_v || !k1_w || !k1_p ||
259 3994 : !k2_u || !k2_v || !k2_w || !k2_p ||
260 3994 : !u0 || !v0 || !w0 || !p0 ||
261 3994 : (needs_T_ws && !T_energy_ws)) {
262 0 : cfd_free(k1_u); cfd_free(k1_v); cfd_free(k1_w); cfd_free(k1_p);
263 0 : cfd_free(k2_u); cfd_free(k2_v); cfd_free(k2_w); cfd_free(k2_p);
264 0 : cfd_free(u0); cfd_free(v0); cfd_free(w0); cfd_free(p0);
265 0 : cfd_free(T_energy_ws);
266 0 : return CFD_ERROR_NOMEM;
267 : }
268 :
269 3994 : double dt = params->dt;
270 3994 : cfd_status_t status = CFD_SUCCESS;
271 :
272 8026 : for (int iter = 0; iter < params->max_iter; iter++) {
273 : /* Save Q^n */
274 4032 : memcpy(u0, field->u, bytes);
275 4032 : memcpy(v0, field->v, bytes);
276 4032 : memcpy(w0, field->w, bytes);
277 4032 : memcpy(p0, field->p, bytes);
278 :
279 : /* ---- Stage 1: k1 = RHS(Q^n) ---- */
280 4032 : memset(k1_u, 0, bytes);
281 4032 : memset(k1_v, 0, bytes);
282 4032 : memset(k1_w, 0, bytes);
283 4032 : memset(k1_p, 0, bytes);
284 :
285 4032 : compute_rhs(field->u, field->v, field->w, field->p, field->rho, field->T,
286 : k1_u, k1_v, k1_w, k1_p,
287 : grid, params, nx, ny, nz,
288 : stride_z, k_start, k_end, inv_2dz, inv_dz2,
289 : iter, dt);
290 :
291 : /* ---- Intermediate: field = Q^n + dt * k1 ---- */
292 44146024 : for (size_t n = 0; n < total; n++) {
293 44141992 : field->u[n] = u0[n] + dt * k1_u[n];
294 44141992 : field->v[n] = v0[n] + dt * k1_v[n];
295 44141992 : field->w[n] = w0[n] + dt * k1_w[n];
296 44141992 : field->p[n] = p0[n] + dt * k1_p[n];
297 :
298 44141992 : field->u[n] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->u[n]));
299 44141992 : field->v[n] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->v[n]));
300 44141992 : field->w[n] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->w[n]));
301 : }
302 :
303 : /* NOTE: Do NOT apply BCs between RK stages. The ghost cells carry
304 : * zero-derivative evolution (k1[ghost]=0), which is consistent with
305 : * the semi-discrete ODE system. Applying BCs here would modify the
306 : * intermediate state outside the ODE trajectory and reduce RK2 to
307 : * first-order temporal accuracy. */
308 :
309 : /* ---- Stage 2: k2 = RHS(Q_pred) ---- */
310 4032 : memset(k2_u, 0, bytes);
311 4032 : memset(k2_v, 0, bytes);
312 4032 : memset(k2_w, 0, bytes);
313 4032 : memset(k2_p, 0, bytes);
314 :
315 4032 : compute_rhs(field->u, field->v, field->w, field->p, field->rho, field->T,
316 : k2_u, k2_v, k2_w, k2_p,
317 : grid, params, nx, ny, nz,
318 : stride_z, k_start, k_end, inv_2dz, inv_dz2,
319 : iter, dt);
320 :
321 : /* ---- Final update: Q^{n+1} = Q^n + (dt/2)*(k1 + k2) ---- */
322 4032 : double half_dt = 0.5 * dt;
323 44146024 : for (size_t n = 0; n < total; n++) {
324 44141992 : field->u[n] = u0[n] + half_dt * (k1_u[n] + k2_u[n]);
325 44141992 : field->v[n] = v0[n] + half_dt * (k1_v[n] + k2_v[n]);
326 44141992 : field->w[n] = w0[n] + half_dt * (k1_w[n] + k2_w[n]);
327 44141992 : field->p[n] = p0[n] + half_dt * (k1_p[n] + k2_p[n]);
328 :
329 44141992 : field->u[n] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->u[n]));
330 44141992 : field->v[n] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->v[n]));
331 44141992 : field->w[n] = fmax(-MAX_VELOCITY_LIMIT, fmin(MAX_VELOCITY_LIMIT, field->w[n]));
332 : }
333 :
334 : /* Energy equation: advance temperature after RK2 velocity update */
335 : {
336 4032 : cfd_status_t energy_status = energy_step_explicit_with_workspace(
337 : field, grid, params, dt, iter * dt, T_energy_ws, total);
338 4032 : if (energy_status != CFD_SUCCESS) {
339 0 : status = energy_status;
340 0 : goto cleanup;
341 : }
342 : }
343 :
344 : /* Apply BCs to final state only (after the full RK2 step).
345 : * This updates ghost cells for the next step's k1 evaluation.
346 : * Then apply configured thermal BCs (overwrites periodic T values). */
347 4032 : apply_boundary_conditions(field, grid);
348 4032 : status = energy_apply_thermal_bcs(field, params);
349 4032 : if (status != CFD_SUCCESS) {
350 0 : goto cleanup;
351 : }
352 :
353 : /* NaN / Inf check */
354 44146024 : for (size_t n = 0; n < total; n++) {
355 44141992 : if (!isfinite(field->u[n]) || !isfinite(field->v[n]) ||
356 44141992 : !isfinite(field->w[n]) || !isfinite(field->p[n])) {
357 0 : status = CFD_ERROR_DIVERGED;
358 0 : goto cleanup;
359 : }
360 : }
361 : }
362 :
363 3994 : cleanup:
364 3994 : cfd_free(k1_u); cfd_free(k1_v); cfd_free(k1_w); cfd_free(k1_p);
365 3994 : cfd_free(k2_u); cfd_free(k2_v); cfd_free(k2_w); cfd_free(k2_p);
366 3994 : cfd_free(u0); cfd_free(v0); cfd_free(w0); cfd_free(p0);
367 3994 : cfd_free(T_energy_ws);
368 :
369 3994 : return status;
370 : }
|
