Line data Source code
1 : /**
2 : * @file energy_solver_omp.c
3 : * @brief OpenMP-parallelized energy equation solver (advection-diffusion step)
4 : *
5 : * Solves: dT/dt + u*nabla(T) = alpha * nabla^2(T) + Q
6 : * using explicit Euler time integration and central differences.
7 : *
8 : * Same numerics as the scalar reference (energy/cpu/energy_solver.c); the
9 : * interior stencil loop is parallelized over j. Read/write separation
10 : * (reads field->T, writes T_new) makes the loop race-free. Buoyancy and
11 : * thermal BC application are shared with the scalar backend.
12 : *
13 : * Branch-free 3D: when nz==1, stride_z=0 and inv_2dz/inv_dz2=0.0 cause all
14 : * z-terms to vanish, producing identical results to a 2D code path.
15 : */
16 :
17 : #include "../energy_solver_internal.h"
18 :
19 : #include "cfd/core/indexing.h"
20 : #include "cfd/core/memory.h"
21 :
22 : #include <math.h>
23 : #include <omp.h>
24 : #include <string.h>
25 :
26 606 : cfd_status_t energy_step_explicit_omp_with_workspace(
27 : flow_field* field, const grid* grid,
28 : const ns_solver_params_t* params,
29 : double dt, double time,
30 : double* T_workspace, size_t workspace_size) {
31 606 : if (!field || !grid || !params) {
32 0 : cfd_set_error(CFD_ERROR_INVALID,
33 : "energy_solver_omp: field, grid, and params must be non-NULL");
34 0 : return CFD_ERROR_INVALID;
35 : }
36 606 : if (!field->T) {
37 0 : cfd_set_error(CFD_ERROR_INVALID,
38 : "energy_solver_omp: missing temperature field");
39 0 : return CFD_ERROR_INVALID;
40 : }
41 :
42 : /* Skip when energy equation is disabled */
43 606 : if (params->alpha <= 0.0) {
44 : return CFD_SUCCESS;
45 : }
46 :
47 60 : size_t nx = field->nx;
48 60 : size_t ny = field->ny;
49 60 : size_t nz = field->nz;
50 :
51 60 : if (!grid->dx || !grid->dy || nx < 3 || ny < 3) {
52 0 : cfd_set_error(CFD_ERROR_INVALID,
53 : "energy_solver_omp: grid too small or missing dx/dy");
54 0 : return CFD_ERROR_INVALID;
55 : }
56 60 : size_t plane = nx * ny;
57 60 : size_t total = plane * nz;
58 60 : double alpha = params->alpha;
59 :
60 : /* Validate uniform spacing (central-difference stencil assumes it) */
61 60 : const double dx0 = grid->dx[0];
62 60 : const double dy0 = grid->dy[0];
63 : {
64 60 : const double tol_x = 1e-12 * fmax(1.0, fabs(dx0));
65 480 : for (size_t i = 1; i < nx - 1; i++) {
66 420 : if (fabs(grid->dx[i] - dx0) > tol_x) {
67 0 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
68 : "energy_solver_omp: non-uniform dx not supported");
69 0 : return CFD_ERROR_UNSUPPORTED;
70 : }
71 : }
72 60 : const double tol_y = 1e-12 * fmax(1.0, fabs(dy0));
73 480 : for (size_t j = 1; j < ny - 1; j++) {
74 420 : if (fabs(grid->dy[j] - dy0) > tol_y) {
75 0 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
76 : "energy_solver_omp: non-uniform dy not supported");
77 0 : return CFD_ERROR_UNSUPPORTED;
78 : }
79 : }
80 : }
81 60 : if (nz > 1) {
82 0 : if (!grid->dz) {
83 0 : cfd_set_error(CFD_ERROR_INVALID,
84 : "energy_solver_omp: missing dz for 3D energy solve");
85 0 : return CFD_ERROR_INVALID;
86 : }
87 0 : const double dz0 = grid->dz[0];
88 0 : const double tol_z = 1e-12 * fmax(1.0, fabs(dz0));
89 0 : for (size_t k = 1; k < nz - 1; k++) {
90 0 : if (fabs(grid->dz[k] - dz0) > tol_z) {
91 0 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
92 : "energy_solver_omp: non-uniform dz not supported");
93 0 : return CFD_ERROR_UNSUPPORTED;
94 : }
95 : }
96 : }
97 :
98 : /* Precomputed constants for uniform-grid stencil */
99 60 : double inv_2dx = 1.0 / (2.0 * dx0);
100 60 : double inv_2dy = 1.0 / (2.0 * dy0);
101 60 : double inv_dx2 = 1.0 / (dx0 * dx0);
102 60 : double inv_dy2 = 1.0 / (dy0 * dy0);
103 :
104 : /* Branch-free 3D constants */
105 60 : size_t stride_z = (nz > 1) ? plane : 0;
106 60 : size_t k_start = (nz > 1) ? 1 : 0;
107 60 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
108 60 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * grid->dz[0]) : 0.0;
109 0 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (grid->dz[0] * grid->dz[0]) : 0.0;
110 :
111 : /* Use caller's workspace or allocate internally */
112 60 : int owns_buffer = 0;
113 60 : double* T_new;
114 60 : if (T_workspace && workspace_size >= total) {
115 : T_new = T_workspace;
116 : } else {
117 0 : T_new = (double*)cfd_calloc(total, sizeof(double));
118 0 : if (!T_new) {
119 : return CFD_ERROR_NOMEM;
120 : }
121 : owns_buffer = 1;
122 : }
123 60 : memcpy(T_new, field->T, total * sizeof(double));
124 :
125 120 : for (size_t k = k_start; k < k_end; k++) {
126 60 : size_t k_offset = k * stride_z;
127 60 : int j;
128 60 : #pragma omp parallel for schedule(static)
129 : for (j = 1; j < (int)ny - 1; j++) {
130 : for (int i = 1; i < (int)nx - 1; i++) {
131 : size_t idx = k_offset + IDX_2D((size_t)i, (size_t)j, nx);
132 :
133 : double T_c = field->T[idx];
134 : double u_c = field->u[idx];
135 : double v_c = field->v[idx];
136 : double w_c = field->w[idx];
137 :
138 : /* Advection: u * dT/dx + v * dT/dy + w * dT/dz */
139 : double dT_dx = (field->T[idx + 1] - field->T[idx - 1]) * inv_2dx;
140 : double dT_dy = (field->T[idx + nx] - field->T[idx - nx]) * inv_2dy;
141 : double dT_dz = (field->T[idx + stride_z] - field->T[idx - stride_z]) * inv_2dz;
142 :
143 : double advection = u_c * dT_dx + v_c * dT_dy + w_c * dT_dz;
144 :
145 : /* Diffusion: alpha * (d2T/dx2 + d2T/dy2 + d2T/dz2) */
146 : double d2T_dx2 = (field->T[idx + 1] - 2.0 * T_c + field->T[idx - 1]) * inv_dx2;
147 : double d2T_dy2 = (field->T[idx + nx] - 2.0 * T_c + field->T[idx - nx]) * inv_dy2;
148 : double d2T_dz2 = (field->T[idx + stride_z] - 2.0 * T_c +
149 : field->T[idx - stride_z]) * inv_dz2;
150 :
151 : double diffusion = alpha * (d2T_dx2 + d2T_dy2 + d2T_dz2);
152 :
153 : /* Heat source term */
154 : double Q = 0.0;
155 : if (params->heat_source_func) {
156 : double x = grid->x[i];
157 : double y = grid->y[j];
158 : double z = (nz > 1 && grid->z) ? grid->z[k] : 0.0;
159 : Q = params->heat_source_func(x, y, z, time,
160 : params->heat_source_context);
161 : }
162 :
163 : /* Explicit Euler update */
164 : double dT = dt * (-advection + diffusion + Q);
165 : T_new[idx] = T_c + dT;
166 : }
167 : }
168 : }
169 :
170 : /* Check for NaN/Inf */
171 60 : int has_nan = 0;
172 60 : ptrdiff_t total_int = (ptrdiff_t)total;
173 60 : ptrdiff_t n;
174 60 : #pragma omp parallel for reduction(| : has_nan) schedule(static)
175 : for (n = 0; n < total_int; n++) {
176 : if (!isfinite(T_new[n])) {
177 : has_nan = 1;
178 : }
179 : }
180 60 : if (has_nan) {
181 0 : cfd_set_error(CFD_ERROR_DIVERGED,
182 : "NaN/Inf detected in energy_step_explicit_omp");
183 0 : if (owns_buffer) cfd_free(T_new);
184 0 : return CFD_ERROR_DIVERGED;
185 : }
186 :
187 60 : memcpy(field->T, T_new, total * sizeof(double));
188 60 : if (owns_buffer) cfd_free(T_new);
189 :
190 : return CFD_SUCCESS;
191 : }
|
