feat: Add per-face thermal BC config to solver params

ROADMAP.md

@@ -26,7 +26,7 @@ This document outlines the development roadmap for achieving a commercial-grade,
 
 ### Backend Coverage Summary
 
-Each algorithm should have scalar (CPU) + SIMD + OMP variants. Track gaps here.
+Each algorithm should have scalar (CPU) + SIMD + OMP + GPU variants. Track gaps here.
 
 | Category            | Algorithm      | CPU  | AVX2     | NEON     | OMP      | GPU  |
 | ------------------- | -------------- | ---- | -------- | -------- | -------- | ---- |
@@ -46,6 +46,7 @@ Each algorithm should have scalar (CPU) + SIMD + OMP variants. Track gaps here.
 - [ ] No turbulence models
 - [ ] Limited linear solvers (no multigrid)
 - [ ] No restart/checkpoint capability
+- [ ] GPU backends limited to projection solver only — missing: Explicit Euler GPU, RK2 GPU, and all modular linear solvers (Jacobi, SOR, Red-Black SOR, CG/PCG, BiCGSTAB)
 
 ### Known Issues
 

lib/include/cfd/solvers/energy_solver.h

@@ -61,6 +61,21 @@ CFD_LIBRARY_EXPORT void energy_compute_buoyancy(double T_local, const ns_solver_
                                                  double* source_u, double* source_v,
                                                  double* source_w);
 
+/**
+ * Apply thermal boundary conditions to the temperature field.
+ *
+ * For each face configured as DIRICHLET, sets T to the specified value.
+ * For each face configured as NEUMANN, copies from adjacent interior cell.
+ * For PERIODIC faces (default), copies from the opposite interior cell.
+ *
+ * Only active when params->alpha > 0. Returns immediately otherwise.
+ *
+ * @param field  Flow field (T is updated in-place)
+ * @param params Solver parameters containing thermal_bc config
+ */
+CFD_LIBRARY_EXPORT void energy_apply_thermal_bcs(flow_field* field,
+                                                  const ns_solver_params_t* params);
+
 #ifdef __cplusplus
 }
 #endif

lib/include/cfd/solvers/navier_stokes_solver.h

@@ -19,6 +19,7 @@
 
 #include "cfd/cfd_export.h"
 
+#include "cfd/boundary/boundary_conditions.h"
 #include "cfd/core/cfd_status.h"
 #include "cfd/core/grid.h"
 #include <stddef.h>
@@ -90,6 +91,26 @@ typedef void (*ns_source_func_t)(double x, double y, double z, double t,
 typedef double (*ns_heat_source_func_t)(double x, double y, double z, double t,
                                          void* context);
 
+/**
+ * Per-face thermal boundary condition configuration.
+ *
+ * Each face can independently be PERIODIC (default), NEUMANN (zero-gradient /
+ * adiabatic), or DIRICHLET (fixed temperature). When a face is DIRICHLET, its
+ * value is read from the corresponding field of `dirichlet_values`.
+ *
+ * Zero-initialization produces all-PERIODIC (BC_TYPE_PERIODIC == 0),
+ * preserving backward compatibility.
+ */
+typedef struct {
+    bc_type_t left;    /**< BC type for x=0 face */
+    bc_type_t right;   /**< BC type for x=Lx face */
+    bc_type_t bottom;  /**< BC type for y=0 face */
+    bc_type_t top;     /**< BC type for y=Ly face */
+    bc_type_t front;   /**< BC type for z=Lz face (3D only) */
+    bc_type_t back;    /**< BC type for z=0 face (3D only) */
+    bc_dirichlet_values_t dirichlet_values;  /**< Fixed values for Dirichlet faces */
+} ns_thermal_bc_config_t;
+
 /**
  * Navier-Stokes solver parameters
  */
@@ -125,6 +146,9 @@ typedef struct {
     /* Custom heat source callback (NULL = no heat source) */
     ns_heat_source_func_t heat_source_func;  /**< Heat source function pointer */
     void* heat_source_context;                /**< User context for heat_source_func */
+
+    /* Thermal boundary conditions (zero-initialized = all PERIODIC = no change) */
+    ns_thermal_bc_config_t thermal_bc;
 } ns_solver_params_t;
 
 

lib/src/solvers/energy/cpu/energy_solver.c

@@ -194,3 +194,111 @@ void energy_compute_buoyancy(double T_local, const ns_solver_params_t* params,
     *source_v += -params->beta * dT * params->gravity[1];
     *source_w += -params->beta * dT * params->gravity[2];
 }
+
+void energy_apply_thermal_bcs(flow_field* field,
+                               const ns_solver_params_t* params) {
+    if (!field || !params || !field->T) return;
+    if (params->alpha <= 0.0) return;
+
+    const ns_thermal_bc_config_t* tbc = &params->thermal_bc;
+
+    size_t nx = field->nx;
+    size_t ny = field->ny;
+    size_t nz = field->nz;
+    size_t plane = nx * ny;
+
+    /* Guard: Neumann needs >= 2 cells, Periodic needs >= 3 cells */
+    if ((tbc->left == BC_TYPE_NEUMANN || tbc->right == BC_TYPE_NEUMANN) && nx < 2) return;
+    if ((tbc->bottom == BC_TYPE_NEUMANN || tbc->top == BC_TYPE_NEUMANN) && ny < 2) return;
+    if ((tbc->left == BC_TYPE_PERIODIC || tbc->right == BC_TYPE_PERIODIC) && nx < 3) return;
+    if ((tbc->bottom == BC_TYPE_PERIODIC || tbc->top == BC_TYPE_PERIODIC) && ny < 3) return;
+    if (nz > 1 && (tbc->back == BC_TYPE_NEUMANN || tbc->front == BC_TYPE_NEUMANN) && nz < 2) return;
+    if (nz > 1 && (tbc->back == BC_TYPE_PERIODIC || tbc->front == BC_TYPE_PERIODIC) && nz < 3) return;
+
+    /* Left face (i=0) */
+    for (size_t k = 0; k < nz; k++) {
+        size_t base = k * plane;
+        for (size_t j = 0; j < ny; j++) {
+            size_t idx = base + j * nx;
+            if (tbc->left == BC_TYPE_DIRICHLET)
+                field->T[idx] = tbc->dirichlet_values.left;
+            else if (tbc->left == BC_TYPE_NEUMANN)
+                field->T[idx] = field->T[idx + 1];
+            else if (tbc->left == BC_TYPE_PERIODIC)
+                field->T[idx] = field->T[base + j * nx + (nx - 2)];
+        }
+    }
+
+    /* Right face (i=nx-1) */
+    for (size_t k = 0; k < nz; k++) {
+        size_t base = k * plane;
+        for (size_t j = 0; j < ny; j++) {
+            size_t idx = base + j * nx + (nx - 1);
+            if (tbc->right == BC_TYPE_DIRICHLET)
+                field->T[idx] = tbc->dirichlet_values.right;
+            else if (tbc->right == BC_TYPE_NEUMANN)
+                field->T[idx] = field->T[idx - 1];
+            else if (tbc->right == BC_TYPE_PERIODIC)
+                field->T[idx] = field->T[base + j * nx + 1];
+        }
+    }
+
+    /* Bottom face (j=0) — runs after left/right, overwrites shared corners */
+    for (size_t k = 0; k < nz; k++) {
+        size_t base = k * plane;
+        for (size_t i = 0; i < nx; i++) {
+            size_t idx = base + i;
+            if (tbc->bottom == BC_TYPE_DIRICHLET)
+                field->T[idx] = tbc->dirichlet_values.bottom;
+            else if (tbc->bottom == BC_TYPE_NEUMANN)
+                field->T[idx] = field->T[idx + nx];
+            else if (tbc->bottom == BC_TYPE_PERIODIC)
+                field->T[idx] = field->T[base + (ny - 2) * nx + i];
+        }
+    }
+
+    /* Top face (j=ny-1) */
+    for (size_t k = 0; k < nz; k++) {
+        size_t base = k * plane;
+        for (size_t i = 0; i < nx; i++) {
+            size_t idx = base + (ny - 1) * nx + i;
+            if (tbc->top == BC_TYPE_DIRICHLET)
+                field->T[idx] = tbc->dirichlet_values.top;
+            else if (tbc->top == BC_TYPE_NEUMANN)
+                field->T[idx] = field->T[idx - nx];
+            else if (tbc->top == BC_TYPE_PERIODIC)
+                field->T[idx] = field->T[base + nx + i];
+        }
+    }
+
+    /* Back face (k=0) — only when nz > 1 */
+    if (nz > 1) {
+        for (size_t j = 0; j < ny; j++) {
+            for (size_t i = 0; i < nx; i++) {
+                size_t idx = j * nx + i;
+                if (tbc->back == BC_TYPE_DIRICHLET)
+                    field->T[idx] = tbc->dirichlet_values.back;
+                else if (tbc->back == BC_TYPE_NEUMANN)
+                    field->T[idx] = field->T[plane + idx];
+                else if (tbc->back == BC_TYPE_PERIODIC)
+                    field->T[idx] = field->T[(nz - 2) * plane + idx];
+            }
+        }
+    }
+
+    /* Front face (k=nz-1) — only when nz > 1 */
+    if (nz > 1) {
+        size_t front_base = (nz - 1) * plane;
+        for (size_t j = 0; j < ny; j++) {
+            for (size_t i = 0; i < nx; i++) {
+                size_t off = j * nx + i;
+                if (tbc->front == BC_TYPE_DIRICHLET)
+                    field->T[front_base + off] = tbc->dirichlet_values.front;
+                else if (tbc->front == BC_TYPE_NEUMANN)
+                    field->T[front_base + off] = field->T[(nz - 2) * plane + off];
+                else if (tbc->front == BC_TYPE_PERIODIC)
+                    field->T[front_base + off] = field->T[plane + off];
+            }
+        }
+    }
+}

lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c

@@ -72,7 +72,8 @@ ns_solver_params_t ns_solver_params_default(void) {
                             .T_ref = 0.0,
                             .gravity = {0.0, 0.0, 0.0},
                             .heat_source_func = NULL,
-                            .heat_source_context = NULL};
+                            .heat_source_context = NULL,
+                            .thermal_bc = {0}};
     return params;
 }
 flow_field* flow_field_create(size_t nx, size_t ny, size_t nz) {
@@ -542,21 +543,14 @@ cfd_status_t explicit_euler_impl(flow_field* field, const grid* grid, const ns_s
             }
         }
 
-        /* Save caller BCs, apply periodic BCs, restore caller BCs.
-         * Use _3d helper for 6-face copy when nz > 1.
-         * Preserve temperature so caller-specified thermal BCs are not
-         * overwritten by the generic periodic BC application. */
+        /* Save caller velocity BCs, apply periodic BCs, restore velocity BCs.
+         * Then apply configured thermal BCs (overwrites periodic T values). */
         copy_boundary_velocities_3d(u_new, v_new, w_new,
                                     field->u, field->v, field->w, nx, ny, nz);
-        if (field->T && T_energy_ws) {
-            memcpy(T_energy_ws, field->T, bytes);
-        }
         apply_boundary_conditions(field, grid);
         copy_boundary_velocities_3d(field->u, field->v, field->w,
                                     u_new, v_new, w_new, nx, ny, nz);
-        if (field->T && T_energy_ws) {
-            memcpy(field->T, T_energy_ws, bytes);
-        }
+        energy_apply_thermal_bcs(field, params);
 
         /* NaN/Inf check */
         int has_nan = 0;

lib/src/solvers/navier_stokes/cpu/solver_projection.c

@@ -264,6 +264,9 @@ cfd_status_t solve_projection_method(flow_field* field, const grid* grid,
             }
         }
 
+        /* Apply configured thermal BCs to temperature field */
+        energy_apply_thermal_bcs(field, params);
+
         /* Restore caller-set boundary conditions */
         copy_boundary_velocities_3d(field->u, field->v, field->w,
                                     u_star, v_star, w_star, nx, ny, nz);

lib/src/solvers/navier_stokes/cpu/solver_rk2.c

@@ -343,15 +343,9 @@ cfd_status_t rk2_impl(flow_field* field, const grid* grid,
 
         /* Apply BCs to final state only (after the full RK2 step).
          * This updates ghost cells for the next step's k1 evaluation.
-         * Preserve temperature so caller-specified thermal BCs are not
-         * overwritten by the generic periodic BC application. */
-        if (field->T && T_energy_ws) {
-            memcpy(T_energy_ws, field->T, bytes);
-            apply_boundary_conditions(field, grid);
-            memcpy(field->T, T_energy_ws, bytes);
-        } else {
-            apply_boundary_conditions(field, grid);
-        }
+         * Then apply configured thermal BCs (overwrites periodic T values). */
+        apply_boundary_conditions(field, grid);
+        energy_apply_thermal_bcs(field, params);
 
         /* NaN / Inf check */
         for (size_t n = 0; n < total; n++) {