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