Line data Source code
1 : /**
2 : * @file energy_solver.c
3 : * @brief Scalar CPU implementation of the energy equation solver
4 : *
5 : * Solves: dT/dt + u*nabla(T) = alpha * nabla^2(T) + Q
6 : * using explicit Euler time integration and central differences.
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 identical results to a 2D code path.
10 : */
11 :
12 : #include "cfd/solvers/energy_solver.h"
13 : #include "../energy_solver_internal.h"
14 :
15 : #include "cfd/core/indexing.h"
16 : #include "cfd/core/memory.h"
17 :
18 : #include <math.h>
19 : #include <string.h>
20 :
21 16602 : cfd_status_t energy_step_explicit_with_workspace(
22 : flow_field* field, const grid* grid,
23 : const ns_solver_params_t* params,
24 : double dt, double time,
25 : double* T_workspace, size_t workspace_size) {
26 16602 : if (!field || !grid || !params) {
27 0 : cfd_set_error(CFD_ERROR_INVALID,
28 : "energy_solver: field, grid, and params must be non-NULL");
29 0 : return CFD_ERROR_INVALID;
30 : }
31 16602 : if (!field->T) {
32 0 : cfd_set_error(CFD_ERROR_INVALID,
33 : "energy_solver: missing temperature field");
34 0 : return CFD_ERROR_INVALID;
35 : }
36 :
37 : /* Skip when energy equation is disabled */
38 16602 : if (params->alpha <= 0.0) {
39 : return CFD_SUCCESS;
40 : }
41 :
42 2199 : size_t nx = field->nx;
43 2199 : size_t ny = field->ny;
44 2199 : size_t nz = field->nz;
45 :
46 2199 : if (!grid->dx || !grid->dy || nx < 3 || ny < 3) {
47 0 : cfd_set_error(CFD_ERROR_INVALID,
48 : "energy_solver: grid too small or missing dx/dy");
49 0 : return CFD_ERROR_INVALID;
50 : }
51 2199 : size_t plane = nx * ny;
52 2199 : size_t total = plane * nz;
53 2199 : double alpha = params->alpha;
54 :
55 : /* Validate uniform spacing (central-difference stencil assumes it) */
56 2199 : const double dx0 = grid->dx[0];
57 2199 : const double dy0 = grid->dy[0];
58 : {
59 2199 : const double tol_x = 1e-12 * fmax(1.0, fabs(dx0));
60 82120 : for (size_t i = 1; i < nx - 1; i++) {
61 79921 : if (fabs(grid->dx[i] - dx0) > tol_x) {
62 0 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
63 : "energy_solver: non-uniform dx not supported");
64 0 : return CFD_ERROR_UNSUPPORTED;
65 : }
66 : }
67 2199 : const double tol_y = 1e-12 * fmax(1.0, fabs(dy0));
68 80320 : for (size_t j = 1; j < ny - 1; j++) {
69 78121 : if (fabs(grid->dy[j] - dy0) > tol_y) {
70 0 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
71 : "energy_solver: non-uniform dy not supported");
72 0 : return CFD_ERROR_UNSUPPORTED;
73 : }
74 : }
75 : }
76 2199 : if (nz > 1) {
77 0 : if (!grid->dz) {
78 0 : cfd_set_error(CFD_ERROR_INVALID,
79 : "energy_solver: missing dz for 3D energy solve");
80 0 : return CFD_ERROR_INVALID;
81 : }
82 0 : const double dz0 = grid->dz[0];
83 0 : const double tol_z = 1e-12 * fmax(1.0, fabs(dz0));
84 0 : for (size_t k = 1; k < nz - 1; k++) {
85 0 : if (fabs(grid->dz[k] - dz0) > tol_z) {
86 0 : cfd_set_error(CFD_ERROR_UNSUPPORTED,
87 : "energy_solver: non-uniform dz not supported");
88 0 : return CFD_ERROR_UNSUPPORTED;
89 : }
90 : }
91 : }
92 :
93 : /* Precomputed constants for uniform-grid stencil */
94 2199 : double inv_2dx = 1.0 / (2.0 * dx0);
95 2199 : double inv_2dy = 1.0 / (2.0 * dy0);
96 2199 : double inv_dx2 = 1.0 / (dx0 * dx0);
97 2199 : double inv_dy2 = 1.0 / (dy0 * dy0);
98 :
99 : /* Branch-free 3D constants */
100 2199 : size_t stride_z = (nz > 1) ? plane : 0;
101 2199 : size_t k_start = (nz > 1) ? 1 : 0;
102 2199 : size_t k_end = (nz > 1) ? (nz - 1) : 1;
103 2199 : double inv_2dz = (nz > 1 && grid->dz) ? 1.0 / (2.0 * grid->dz[0]) : 0.0;
104 0 : double inv_dz2 = (nz > 1 && grid->dz) ? 1.0 / (grid->dz[0] * grid->dz[0]) : 0.0;
105 :
106 : /* Use caller's workspace or allocate internally */
107 2199 : int owns_buffer = 0;
108 2199 : double* T_new;
109 2199 : if (T_workspace && workspace_size >= total) {
110 : T_new = T_workspace;
111 : } else {
112 200 : T_new = (double*)cfd_calloc(total, sizeof(double));
113 200 : if (!T_new) {
114 : return CFD_ERROR_NOMEM;
115 : }
116 : owns_buffer = 1;
117 : }
118 2199 : memcpy(T_new, field->T, total * sizeof(double));
119 :
120 4398 : for (size_t k = k_start; k < k_end; k++) {
121 80320 : for (size_t j = 1; j < ny - 1; j++) {
122 3055160 : for (size_t i = 1; i < nx - 1; i++) {
123 2977039 : size_t idx = k * stride_z + IDX_2D(i, j, nx);
124 :
125 2977039 : double T_c = field->T[idx];
126 2977039 : double u_c = field->u[idx];
127 2977039 : double v_c = field->v[idx];
128 2977039 : double w_c = field->w[idx];
129 :
130 : /* Advection: u * dT/dx + v * dT/dy + w * dT/dz */
131 2977039 : double dT_dx = (field->T[idx + 1] - field->T[idx - 1]) * inv_2dx;
132 2977039 : double dT_dy = (field->T[idx + nx] - field->T[idx - nx]) * inv_2dy;
133 2977039 : double dT_dz = (field->T[idx + stride_z] - field->T[idx - stride_z]) * inv_2dz;
134 :
135 2977039 : double advection = u_c * dT_dx + v_c * dT_dy + w_c * dT_dz;
136 :
137 : /* Diffusion: alpha * (d2T/dx2 + d2T/dy2 + d2T/dz2) */
138 2977039 : double d2T_dx2 = (field->T[idx + 1] - 2.0 * T_c + field->T[idx - 1]) * inv_dx2;
139 2977039 : double d2T_dy2 = (field->T[idx + nx] - 2.0 * T_c + field->T[idx - nx]) * inv_dy2;
140 2977039 : double d2T_dz2 = (field->T[idx + stride_z] - 2.0 * T_c +
141 : field->T[idx - stride_z]) * inv_dz2;
142 :
143 2977039 : double diffusion = alpha * (d2T_dx2 + d2T_dy2 + d2T_dz2);
144 :
145 : /* Heat source term */
146 2977039 : double Q = 0.0;
147 2977039 : if (params->heat_source_func) {
148 11250 : double x = grid->x[i];
149 11250 : double y = grid->y[j];
150 11250 : double z = (nz > 1 && grid->z) ? grid->z[k] : 0.0;
151 11250 : Q = params->heat_source_func(x, y, z, time,
152 11250 : params->heat_source_context);
153 : }
154 :
155 : /* Explicit Euler update */
156 2977039 : double dT = dt * (-advection + diffusion + Q);
157 2977039 : T_new[idx] = T_c + dT;
158 : }
159 : }
160 : }
161 :
162 : /* Check for NaN/Inf */
163 3304118 : for (size_t n = 0; n < total; n++) {
164 3301919 : if (!isfinite(T_new[n])) {
165 0 : cfd_set_error(CFD_ERROR_DIVERGED,
166 : "NaN/Inf detected in energy_step_explicit");
167 0 : if (owns_buffer) cfd_free(T_new);
168 0 : return CFD_ERROR_DIVERGED;
169 : }
170 : }
171 :
172 2199 : memcpy(field->T, T_new, total * sizeof(double));
173 2199 : if (owns_buffer) cfd_free(T_new);
174 :
175 : return CFD_SUCCESS;
176 : }
177 :
178 201 : cfd_status_t energy_step_explicit(flow_field* field, const grid* grid,
179 : const ns_solver_params_t* params,
180 : double dt, double time) {
181 201 : return energy_step_explicit_with_workspace(field, grid, params,
182 : dt, time, NULL, 0);
183 : }
184 :
185 141162014 : void energy_compute_buoyancy(double T_local, const ns_solver_params_t* params,
186 : double* source_u, double* source_v,
187 : double* source_w) {
188 141162014 : if (params->beta == 0.0) {
189 : return;
190 : }
191 :
192 2857961 : double dT = T_local - params->T_ref;
193 2857961 : *source_u += -params->beta * dT * params->gravity[0];
194 2857961 : *source_v += -params->beta * dT * params->gravity[1];
195 2857961 : *source_w += -params->beta * dT * params->gravity[2];
196 : }
197 :
198 : /* A face may only request a thermal BC type the energy solver implements. */
199 : static int is_supported_thermal_bc(bc_type_t type) {
200 : return type == BC_TYPE_PERIODIC || type == BC_TYPE_NEUMANN ||
201 : type == BC_TYPE_DIRICHLET;
202 : }
203 :
204 17340 : cfd_status_t energy_apply_thermal_bcs(flow_field* field,
205 : const ns_solver_params_t* params) {
206 17340 : if (!field || !params || !field->T) {
207 2 : cfd_set_error(CFD_ERROR_INVALID,
208 : "energy_apply_thermal_bcs: field, params, and T must be non-NULL");
209 2 : return CFD_ERROR_INVALID;
210 : }
211 : /* alpha <= 0 means the energy equation is disabled: a legitimate no-op. */
212 17338 : if (params->alpha <= 0.0) return CFD_SUCCESS;
213 :
214 2084 : const ns_thermal_bc_config_t* tbc = ¶ms->thermal_bc;
215 :
216 2084 : size_t nx = field->nx;
217 2084 : size_t ny = field->ny;
218 2084 : size_t nz = field->nz;
219 2084 : size_t plane = nx * ny;
220 :
221 : /* Reject unsupported per-face BC types so misconfiguration (e.g. an
222 : * accidental BC_TYPE_NOSLIP/INLET) fails loudly instead of silently
223 : * leaving that face unchanged. Front/back only apply in 3D. */
224 2084 : if (!is_supported_thermal_bc(tbc->left) || !is_supported_thermal_bc(tbc->right) ||
225 2083 : !is_supported_thermal_bc(tbc->bottom) || !is_supported_thermal_bc(tbc->top) ||
226 1 : (nz > 1 && (!is_supported_thermal_bc(tbc->front) ||
227 1 : !is_supported_thermal_bc(tbc->back)))) {
228 1 : cfd_set_error(CFD_ERROR_INVALID,
229 : "energy_apply_thermal_bcs: unsupported thermal BC type on a face "
230 : "(only PERIODIC, NEUMANN, DIRICHLET are valid)");
231 1 : return CFD_ERROR_INVALID;
232 : }
233 :
234 : /* A requested BC must fit the grid: Neumann needs >= 2 cells, Periodic >= 3. */
235 2083 : if (((tbc->left == BC_TYPE_NEUMANN || tbc->right == BC_TYPE_NEUMANN) && nx < 2) ||
236 2083 : ((tbc->bottom == BC_TYPE_NEUMANN || tbc->top == BC_TYPE_NEUMANN) && ny < 2) ||
237 2083 : ((tbc->left == BC_TYPE_PERIODIC || tbc->right == BC_TYPE_PERIODIC) && nx < 3) ||
238 2083 : ((tbc->bottom == BC_TYPE_PERIODIC || tbc->top == BC_TYPE_PERIODIC) && ny < 3) ||
239 2083 : (nz > 1 && (tbc->back == BC_TYPE_NEUMANN || tbc->front == BC_TYPE_NEUMANN) && nz < 2) ||
240 1 : (nz > 1 && (tbc->back == BC_TYPE_PERIODIC || tbc->front == BC_TYPE_PERIODIC) && nz < 3)) {
241 0 : cfd_set_error(CFD_ERROR_INVALID,
242 : "energy_apply_thermal_bcs: grid too small for the requested thermal BC type");
243 0 : return CFD_ERROR_INVALID;
244 : }
245 :
246 : /* Left face (i=0) */
247 4170 : for (size_t k = 0; k < nz; k++) {
248 2087 : size_t base = k * plane;
249 80986 : for (size_t j = 0; j < ny; j++) {
250 78899 : size_t idx = base + j * nx;
251 78899 : if (tbc->left == BC_TYPE_DIRICHLET)
252 78855 : field->T[idx] = tbc->dirichlet_values.left;
253 44 : else if (tbc->left == BC_TYPE_NEUMANN)
254 0 : field->T[idx] = field->T[idx + 1];
255 44 : else if (tbc->left == BC_TYPE_PERIODIC)
256 44 : field->T[idx] = field->T[base + j * nx + (nx - 2)];
257 : }
258 : }
259 :
260 : /* Right face (i=nx-1) */
261 4170 : for (size_t k = 0; k < nz; k++) {
262 2087 : size_t base = k * plane;
263 80986 : for (size_t j = 0; j < ny; j++) {
264 78899 : size_t idx = base + j * nx + (nx - 1);
265 78899 : if (tbc->right == BC_TYPE_DIRICHLET)
266 78850 : field->T[idx] = tbc->dirichlet_values.right;
267 49 : else if (tbc->right == BC_TYPE_NEUMANN)
268 0 : field->T[idx] = field->T[idx - 1];
269 49 : else if (tbc->right == BC_TYPE_PERIODIC)
270 49 : field->T[idx] = field->T[base + j * nx + 1];
271 : }
272 : }
273 :
274 : /* Bottom face (j=0) — runs after left/right, overwrites shared corners */
275 4170 : for (size_t k = 0; k < nz; k++) {
276 2087 : size_t base = k * plane;
277 80986 : for (size_t i = 0; i < nx; i++) {
278 78899 : size_t idx = base + i;
279 78899 : if (tbc->bottom == BC_TYPE_DIRICHLET)
280 35 : field->T[idx] = tbc->dirichlet_values.bottom;
281 78864 : else if (tbc->bottom == BC_TYPE_NEUMANN)
282 78850 : field->T[idx] = field->T[idx + nx];
283 14 : else if (tbc->bottom == BC_TYPE_PERIODIC)
284 14 : field->T[idx] = field->T[base + (ny - 2) * nx + i];
285 : }
286 : }
287 :
288 : /* Top face (j=ny-1) */
289 4170 : for (size_t k = 0; k < nz; k++) {
290 2087 : size_t base = k * plane;
291 80986 : for (size_t i = 0; i < nx; i++) {
292 78899 : size_t idx = base + (ny - 1) * nx + i;
293 78899 : if (tbc->top == BC_TYPE_DIRICHLET)
294 35 : field->T[idx] = tbc->dirichlet_values.top;
295 78864 : else if (tbc->top == BC_TYPE_NEUMANN)
296 78850 : field->T[idx] = field->T[idx - nx];
297 14 : else if (tbc->top == BC_TYPE_PERIODIC)
298 14 : field->T[idx] = field->T[base + nx + i];
299 : }
300 : }
301 :
302 : /* Back face (k=0) — only when nz > 1 */
303 2083 : if (nz > 1) {
304 8 : for (size_t j = 0; j < ny; j++) {
305 56 : for (size_t i = 0; i < nx; i++) {
306 49 : size_t idx = j * nx + i;
307 49 : if (tbc->back == BC_TYPE_DIRICHLET)
308 49 : field->T[idx] = tbc->dirichlet_values.back;
309 0 : else if (tbc->back == BC_TYPE_NEUMANN)
310 0 : field->T[idx] = field->T[plane + idx];
311 0 : else if (tbc->back == BC_TYPE_PERIODIC)
312 0 : field->T[idx] = field->T[(nz - 2) * plane + idx];
313 : }
314 : }
315 : }
316 :
317 : /* Front face (k=nz-1) — only when nz > 1 */
318 1 : if (nz > 1) {
319 1 : size_t front_base = (nz - 1) * plane;
320 8 : for (size_t j = 0; j < ny; j++) {
321 56 : for (size_t i = 0; i < nx; i++) {
322 49 : size_t off = j * nx + i;
323 49 : if (tbc->front == BC_TYPE_DIRICHLET)
324 0 : field->T[front_base + off] = tbc->dirichlet_values.front;
325 49 : else if (tbc->front == BC_TYPE_NEUMANN)
326 49 : field->T[front_base + off] = field->T[(nz - 2) * plane + off];
327 0 : else if (tbc->front == BC_TYPE_PERIODIC)
328 0 : field->T[front_base + off] = field->T[plane + off];
329 : }
330 : }
331 : }
332 :
333 : return CFD_SUCCESS;
334 : }
|
