From ed5d256ca5acb4d73b18caf585241a3a4053748b Mon Sep 17 00:00:00 2001 From: shaia Date: Sun, 5 Apr 2026 00:44:53 +0300 Subject: [PATCH 1/5] feat: Add per-face thermal BC config to solver params Thermal boundary conditions (Dirichlet/Neumann) are now configured via ns_solver_params_t.thermal_bc instead of being applied manually by the caller. The solver applies them automatically after each energy step. - Add ns_thermal_bc_config_t with per-face bc_type_t + Dirichlet values - Add energy_apply_thermal_bcs() for 3D-aware thermal BC application - Remove save/restore T memcpy hack from explicit Euler, projection, RK2 - Update natural convection test to use config-driven thermal BCs - Add unit tests for mixed Dirichlet/Neumann and all-periodic no-op --- ROADMAP.md | 3 +- lib/include/cfd/solvers/energy_solver.h | 16 +++ .../cfd/solvers/navier_stokes_solver.h | 24 +++++ lib/src/solvers/energy/cpu/energy_solver.c | 96 +++++++++++++++++ .../navier_stokes/cpu/solver_explicit_euler.c | 16 +-- .../navier_stokes/cpu/solver_projection.c | 3 + .../solvers/navier_stokes/cpu/solver_rk2.c | 12 +-- tests/solvers/energy/test_energy_solver.c | 102 ++++++++++++++++++ tests/validation/test_natural_convection.c | 38 +++---- 9 files changed, 268 insertions(+), 42 deletions(-) diff --git a/ROADMAP.md b/ROADMAP.md index 8212690..84f4825 100644 --- a/ROADMAP.md +++ b/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 diff --git a/lib/include/cfd/solvers/energy_solver.h b/lib/include/cfd/solvers/energy_solver.h index d632435..02dcf77 100644 --- a/lib/include/cfd/solvers/energy_solver.h +++ b/lib/include/cfd/solvers/energy_solver.h @@ -61,6 +61,22 @@ 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), does nothing (assumes periodic BCs were + * already applied by apply_boundary_conditions). + * + * 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 diff --git a/lib/include/cfd/solvers/navier_stokes_solver.h b/lib/include/cfd/solvers/navier_stokes_solver.h index 9df26c4..1e26504 100644 --- a/lib/include/cfd/solvers/navier_stokes_solver.h +++ b/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 @@ -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; diff --git a/lib/src/solvers/energy/cpu/energy_solver.c b/lib/src/solvers/energy/cpu/energy_solver.c index 154fd94..4bc343e 100644 --- a/lib/src/solvers/energy/cpu/energy_solver.c +++ b/lib/src/solvers/energy/cpu/energy_solver.c @@ -194,3 +194,99 @@ 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 = ¶ms->thermal_bc; + + /* All-periodic means nothing to do */ + if (tbc->left == BC_TYPE_PERIODIC && tbc->right == BC_TYPE_PERIODIC && + tbc->bottom == BC_TYPE_PERIODIC && tbc->top == BC_TYPE_PERIODIC && + tbc->front == BC_TYPE_PERIODIC && tbc->back == BC_TYPE_PERIODIC) { + return; + } + + size_t nx = field->nx; + size_t ny = field->ny; + size_t nz = field->nz; + size_t plane = nx * ny; + + /* 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]; + } + } + + /* 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]; + } + } + + /* 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]; + } + } + + /* 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]; + } + } + + /* 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]; + } + } + } + + /* Front face (k=nz-1) — only when nz > 1 */ + if (nz > 1) { + size_t front_base = (nz - 1) * plane; + size_t interior_base = (nz - 2) * 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[interior_base + off]; + } + } + } +} diff --git a/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c b/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c index 0049068..0561db9 100644 --- a/lib/src/solvers/navier_stokes/cpu/solver_explicit_euler.c +++ b/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; diff --git a/lib/src/solvers/navier_stokes/cpu/solver_projection.c b/lib/src/solvers/navier_stokes/cpu/solver_projection.c index d73fa47..2e3020f 100644 --- a/lib/src/solvers/navier_stokes/cpu/solver_projection.c +++ b/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); diff --git a/lib/src/solvers/navier_stokes/cpu/solver_rk2.c b/lib/src/solvers/navier_stokes/cpu/solver_rk2.c index 58d5efc..dc24ba2 100644 --- a/lib/src/solvers/navier_stokes/cpu/solver_rk2.c +++ b/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++) { diff --git a/tests/solvers/energy/test_energy_solver.c b/tests/solvers/energy/test_energy_solver.c index 418b157..16b2d36 100644 --- a/tests/solvers/energy/test_energy_solver.c +++ b/tests/solvers/energy/test_energy_solver.c @@ -439,6 +439,106 @@ static void test_heat_source(void) { grid_destroy(g); } +/* ============================================================================ + * TEST 7: Thermal BC — Mixed Dirichlet / Neumann + * + * Configure Dirichlet on left/right (fixed T) and Neumann on top/bottom + * (zero-gradient). After calling energy_apply_thermal_bcs, verify: + * - Left/right boundaries have the Dirichlet values + * - Top/bottom boundaries equal the adjacent interior row + * ============================================================================ */ + +#define TBC_NX 11 +#define TBC_NY 11 + +static void test_thermal_bc_dirichlet_neumann(void) { + grid* g = grid_create(TBC_NX, TBC_NY, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + flow_field* field = flow_field_create(TBC_NX, TBC_NY, 1); + TEST_ASSERT_NOT_NULL(field); + + /* Fill T with a recognizable pattern: T = 100 + i + j*10 */ + for (size_t j = 0; j < TBC_NY; j++) { + for (size_t i = 0; i < TBC_NX; i++) { + size_t idx = j * TBC_NX + i; + field->T[idx] = 100.0 + (double)i + (double)j * 10.0; + field->rho[idx] = 1.0; + } + } + + ns_solver_params_t params = ns_solver_params_default(); + params.alpha = 0.01; /* Must be > 0 to activate thermal BCs */ + params.thermal_bc.left = BC_TYPE_DIRICHLET; + params.thermal_bc.right = BC_TYPE_DIRICHLET; + params.thermal_bc.bottom = BC_TYPE_NEUMANN; + params.thermal_bc.top = BC_TYPE_NEUMANN; + params.thermal_bc.dirichlet_values.left = 500.0; + params.thermal_bc.dirichlet_values.right = 200.0; + + energy_apply_thermal_bcs(field, ¶ms); + + /* Verify Dirichlet on left (i=0) and right (i=nx-1) */ + for (size_t j = 0; j < TBC_NY; j++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 500.0, field->T[j * TBC_NX + 0]); + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 200.0, field->T[j * TBC_NX + (TBC_NX - 1)]); + } + + /* Verify Neumann on bottom (j=0): T[0,i] == T[1,i] */ + for (size_t i = 0; i < TBC_NX; i++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, field->T[1 * TBC_NX + i], + field->T[0 * TBC_NX + i]); + } + + /* Verify Neumann on top (j=ny-1): T[ny-1,i] == T[ny-2,i] */ + for (size_t i = 0; i < TBC_NX; i++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, field->T[(TBC_NY - 2) * TBC_NX + i], + field->T[(TBC_NY - 1) * TBC_NX + i]); + } + + flow_field_destroy(field); + grid_destroy(g); +} + +/* ============================================================================ + * TEST 8: Thermal BC — All Periodic (no-op) + * + * Default thermal_bc config (all PERIODIC) should not modify T at all. + * ============================================================================ */ + +static void test_thermal_bc_all_periodic_noop(void) { + grid* g = grid_create(9, 9, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + flow_field* field = flow_field_create(9, 9, 1); + TEST_ASSERT_NOT_NULL(field); + + /* Fill T with distinct values */ + for (size_t n = 0; n < 81; n++) { + field->T[n] = 1000.0 + (double)n; + field->rho[n] = 1.0; + } + + /* Save original T */ + double T_orig[81]; + memcpy(T_orig, field->T, 81 * sizeof(double)); + + ns_solver_params_t params = ns_solver_params_default(); + params.alpha = 0.01; /* Energy enabled, but thermal_bc is all-periodic */ + + energy_apply_thermal_bcs(field, ¶ms); + + /* T should be completely unchanged */ + for (size_t n = 0; n < 81; n++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[n], field->T[n]); + } + + flow_field_destroy(field); + grid_destroy(g); +} + /* ============================================================================ * MAIN * ============================================================================ */ @@ -451,5 +551,7 @@ int main(void) { RUN_TEST(test_buoyancy_source); RUN_TEST(test_backward_compatibility); RUN_TEST(test_heat_source); + RUN_TEST(test_thermal_bc_dirichlet_neumann); + RUN_TEST(test_thermal_bc_all_periodic_noop); return UNITY_END(); } diff --git a/tests/validation/test_natural_convection.c b/tests/validation/test_natural_convection.c index 2ff959e..79b878e 100644 --- a/tests/validation/test_natural_convection.c +++ b/tests/validation/test_natural_convection.c @@ -22,6 +22,7 @@ * 4. No divergence (NaN/Inf) */ +#include "cfd/boundary/boundary_conditions.h" #include "cfd/core/cfd_init.h" #include "cfd/core/cfd_status.h" #include "cfd/core/grid.h" @@ -70,33 +71,19 @@ void setUp(void) { cfd_init(); } void tearDown(void) { cfd_finalize(); } /* ============================================================================ - * Apply thermal boundary conditions for differentially heated cavity + * Apply velocity no-slip boundary conditions for cavity walls + * (Thermal BCs are now configured via params.thermal_bc) * ============================================================================ */ -static void apply_cavity_bcs(flow_field* field, size_t nx, size_t ny) { +static void apply_cavity_velocity_bcs(flow_field* field, size_t nx, size_t ny) { for (size_t j = 0; j < ny; j++) { - /* Left wall: hot */ - field->T[j * nx + 0] = NC_T_HOT; - /* Right wall: cold */ - field->T[j * nx + (nx - 1)] = NC_T_COLD; - - /* No-slip on left and right walls */ field->u[j * nx + 0] = 0.0; field->v[j * nx + 0] = 0.0; field->u[j * nx + (nx - 1)] = 0.0; field->v[j * nx + (nx - 1)] = 0.0; } - for (size_t i = 1; i < nx - 1; i++) { - /* Bottom wall: adiabatic (Neumann dT/dy=0 via copy from interior) */ - field->T[0 * nx + i] = field->T[1 * nx + i]; - /* Top wall: adiabatic */ - field->T[(ny - 1) * nx + i] = field->T[(ny - 2) * nx + i]; - } - for (size_t i = 0; i < nx; i++) { - - /* No-slip on top and bottom walls */ field->u[0 * nx + i] = 0.0; field->v[0 * nx + i] = 0.0; field->u[(ny - 1) * nx + i] = 0.0; @@ -145,6 +132,14 @@ static void test_natural_convection_development(void) { params.source_amplitude_u = 0.0; params.source_amplitude_v = 0.0; + /* Thermal BCs: Dirichlet on left/right walls, Neumann (adiabatic) on top/bottom */ + params.thermal_bc.left = BC_TYPE_DIRICHLET; + params.thermal_bc.right = BC_TYPE_DIRICHLET; + params.thermal_bc.top = BC_TYPE_NEUMANN; + params.thermal_bc.bottom = BC_TYPE_NEUMANN; + params.thermal_bc.dirichlet_values.left = NC_T_HOT; + params.thermal_bc.dirichlet_values.right = NC_T_COLD; + /* Use projection solver */ ns_solver_registry_t* registry = cfd_registry_create(); TEST_ASSERT_NOT_NULL(registry); @@ -159,15 +154,16 @@ static void test_natural_convection_development(void) { /* Run simulation */ ns_solver_stats_t stats = ns_solver_stats_default(); for (int step = 0; step < NC_STEPS; step++) { - /* Apply thermal BCs before each step */ - apply_cavity_bcs(field, NC_NX, NC_NY); + /* Apply velocity no-slip BCs before step */ + apply_cavity_velocity_bcs(field, NC_NX, NC_NY); cfd_status_t status = solver_step(solver, field, g, ¶ms, &stats); TEST_ASSERT_EQUAL_MESSAGE(CFD_SUCCESS, status, "Solver step should succeed"); - /* Re-apply thermal BCs after step */ - apply_cavity_bcs(field, NC_NX, NC_NY); + /* Thermal BCs are applied automatically inside solver_step. + * Reapply velocity no-slip after step. */ + apply_cavity_velocity_bcs(field, NC_NX, NC_NY); } /* ---- Verify results ---- */ From d78ee8e47a4b8a741ec5caddc33ff8d9a3ee57c5 Mon Sep 17 00:00:00 2001 From: shaia Date: Tue, 7 Apr 2026 20:34:44 +0300 Subject: [PATCH 2/5] Add guards and periodic handling to energy_apply_thermal_bcs Neumann BCs read from adjacent interior cells, which is out-of-bounds when nx<2 or ny<2. Add early-return guards for undersized grids. Periodic temperature BCs were previously a no-op, relying on the caller to have already called apply_boundary_conditions(). This left periodic temperature boundaries unenforced in the projection solver. Handle all three BC types (Dirichlet, Neumann, Periodic) directly so the function is self-contained across all solvers. --- lib/include/cfd/solvers/energy_solver.h | 3 +-- lib/src/solvers/energy/cpu/energy_solver.c | 30 +++++++++++++++------- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/lib/include/cfd/solvers/energy_solver.h b/lib/include/cfd/solvers/energy_solver.h index 02dcf77..663d299 100644 --- a/lib/include/cfd/solvers/energy_solver.h +++ b/lib/include/cfd/solvers/energy_solver.h @@ -66,8 +66,7 @@ CFD_LIBRARY_EXPORT void energy_compute_buoyancy(double T_local, const ns_solver_ * * 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), does nothing (assumes periodic BCs were - * already applied by apply_boundary_conditions). + * For PERIODIC faces (default), copies from the opposite interior cell. * * Only active when params->alpha > 0. Returns immediately otherwise. * diff --git a/lib/src/solvers/energy/cpu/energy_solver.c b/lib/src/solvers/energy/cpu/energy_solver.c index 4bc343e..5303746 100644 --- a/lib/src/solvers/energy/cpu/energy_solver.c +++ b/lib/src/solvers/energy/cpu/energy_solver.c @@ -202,18 +202,19 @@ void energy_apply_thermal_bcs(flow_field* field, const ns_thermal_bc_config_t* tbc = ¶ms->thermal_bc; - /* All-periodic means nothing to do */ - if (tbc->left == BC_TYPE_PERIODIC && tbc->right == BC_TYPE_PERIODIC && - tbc->bottom == BC_TYPE_PERIODIC && tbc->top == BC_TYPE_PERIODIC && - tbc->front == BC_TYPE_PERIODIC && tbc->back == BC_TYPE_PERIODIC) { - return; - } - 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; @@ -223,6 +224,8 @@ void energy_apply_thermal_bcs(flow_field* field, 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)]; } } @@ -235,6 +238,8 @@ void energy_apply_thermal_bcs(flow_field* field, 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]; } } @@ -247,6 +252,8 @@ void energy_apply_thermal_bcs(flow_field* field, 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]; } } @@ -259,6 +266,8 @@ void energy_apply_thermal_bcs(flow_field* field, 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]; } } @@ -271,6 +280,8 @@ void energy_apply_thermal_bcs(flow_field* field, 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]; } } } @@ -278,14 +289,15 @@ void energy_apply_thermal_bcs(flow_field* field, /* Front face (k=nz-1) — only when nz > 1 */ if (nz > 1) { size_t front_base = (nz - 1) * plane; - size_t interior_base = (nz - 2) * 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[interior_base + off]; + 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]; } } } From e6cbe69af73f3def0ee3c15acc3d78834120ec88 Mon Sep 17 00:00:00 2001 From: shaia Date: Tue, 7 Apr 2026 20:37:48 +0300 Subject: [PATCH 3/5] Add 3D thermal BC test for front/back and periodic faces Exercises nz>1 code paths: Dirichlet on back face, Neumann on front face, Periodic on left/right, and Dirichlet on bottom/top. Verifies corner overwrite ordering and periodic copy from opposite interior. --- tests/solvers/energy/test_energy_solver.c | 127 ++++++++++++++++++++++ 1 file changed, 127 insertions(+) diff --git a/tests/solvers/energy/test_energy_solver.c b/tests/solvers/energy/test_energy_solver.c index 16b2d36..494e553 100644 --- a/tests/solvers/energy/test_energy_solver.c +++ b/tests/solvers/energy/test_energy_solver.c @@ -539,6 +539,132 @@ static void test_thermal_bc_all_periodic_noop(void) { grid_destroy(g); } +/* ============================================================================ + * TEST 9: 3D Thermal BCs — Front/Back Dirichlet + Neumann + * + * Create a small 3D field (nz=5) and verify: + * - Back face (k=0): Dirichlet sets T to specified value + * - Front face (k=nz-1): Neumann copies from adjacent interior plane + * - Left/right: Periodic copies from opposite interior column + * - Bottom/top: Dirichlet values overwrite corners correctly + * ============================================================================ */ + +#define TBC3D_NX 7 +#define TBC3D_NY 7 +#define TBC3D_NZ 5 +#define TBC3D_PLANE ((size_t)(TBC3D_NX) * TBC3D_NY) + +static void test_thermal_bc_3d_front_back(void) { + grid* g = grid_create(TBC3D_NX, TBC3D_NY, TBC3D_NZ, + 0.0, 1.0, 0.0, 1.0, 0.0, 1.0); + TEST_ASSERT_NOT_NULL(g); + grid_initialize_uniform(g); + + flow_field* field = flow_field_create(TBC3D_NX, TBC3D_NY, TBC3D_NZ); + TEST_ASSERT_NOT_NULL(field); + + size_t total = TBC3D_PLANE * TBC3D_NZ; + + /* Fill T with a recognizable 3D pattern: T = k*1000 + j*100 + i */ + for (size_t k = 0; k < TBC3D_NZ; k++) { + for (size_t j = 0; j < TBC3D_NY; j++) { + for (size_t i = 0; i < TBC3D_NX; i++) { + size_t idx = k * TBC3D_PLANE + j * TBC3D_NX + i; + field->T[idx] = (double)k * 1000.0 + (double)j * 100.0 + (double)i; + field->rho[idx] = 1.0; + } + } + } + + ns_solver_params_t params = ns_solver_params_default(); + params.alpha = 0.01; /* Enable thermal BCs */ + + /* Back = Dirichlet(50), Front = Neumann, Left/Right = Periodic, + * Bottom = Dirichlet(99), Top = Dirichlet(77) */ + params.thermal_bc.back = BC_TYPE_DIRICHLET; + params.thermal_bc.front = BC_TYPE_NEUMANN; + params.thermal_bc.left = BC_TYPE_PERIODIC; + params.thermal_bc.right = BC_TYPE_PERIODIC; + params.thermal_bc.bottom = BC_TYPE_DIRICHLET; + params.thermal_bc.top = BC_TYPE_DIRICHLET; + params.thermal_bc.dirichlet_values.back = 50.0; + params.thermal_bc.dirichlet_values.bottom = 99.0; + params.thermal_bc.dirichlet_values.top = 77.0; + + energy_apply_thermal_bcs(field, ¶ms); + + /* Verify back face (k=0): all cells should be 50.0 (Dirichlet) */ + for (size_t j = 0; j < TBC3D_NY; j++) { + for (size_t i = 0; i < TBC3D_NX; i++) { + size_t idx = j * TBC3D_NX + i; + /* Bottom/top Dirichlet overwrites corners on j=0 and j=ny-1 */ + if (j == 0) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 99.0, field->T[idx]); + } else if (j == TBC3D_NY - 1) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 77.0, field->T[idx]); + } else { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 50.0, field->T[idx]); + } + } + } + + /* Verify front face (k=nz-1): Neumann copies from k=nz-2 */ + size_t front_base = (TBC3D_NZ - 1) * TBC3D_PLANE; + size_t interior_base = (TBC3D_NZ - 2) * TBC3D_PLANE; + for (size_t j = 1; j < TBC3D_NY - 1; j++) { + for (size_t i = 1; i < TBC3D_NX - 1; i++) { + size_t off = j * TBC3D_NX + i; + TEST_ASSERT_DOUBLE_WITHIN(1e-15, + field->T[interior_base + off], + field->T[front_base + off]); + } + } + + /* Verify left face (i=0): Periodic copies from i=nx-2 */ + for (size_t k = 1; k < TBC3D_NZ - 1; k++) { + size_t base = k * TBC3D_PLANE; + for (size_t j = 1; j < TBC3D_NY - 1; j++) { + size_t left_idx = base + j * TBC3D_NX; + size_t src_idx = base + j * TBC3D_NX + (TBC3D_NX - 2); + TEST_ASSERT_DOUBLE_WITHIN(1e-15, + field->T[src_idx], + field->T[left_idx]); + } + } + + /* Verify right face (i=nx-1): Periodic copies from i=1 */ + for (size_t k = 1; k < TBC3D_NZ - 1; k++) { + size_t base = k * TBC3D_PLANE; + for (size_t j = 1; j < TBC3D_NY - 1; j++) { + size_t right_idx = base + j * TBC3D_NX + (TBC3D_NX - 1); + size_t src_idx = base + j * TBC3D_NX + 1; + TEST_ASSERT_DOUBLE_WITHIN(1e-15, + field->T[src_idx], + field->T[right_idx]); + } + } + + /* Verify bottom face (j=0): Dirichlet 99.0 across all k */ + for (size_t k = 0; k < TBC3D_NZ; k++) { + size_t base = k * TBC3D_PLANE; + for (size_t i = 0; i < TBC3D_NX; i++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 99.0, field->T[base + i]); + } + } + + /* Verify top face (j=ny-1): Dirichlet 77.0 across all k */ + for (size_t k = 0; k < TBC3D_NZ; k++) { + size_t base = k * TBC3D_PLANE; + for (size_t i = 0; i < TBC3D_NX; i++) { + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 77.0, + field->T[base + (TBC3D_NY - 1) * TBC3D_NX + i]); + } + } + + flow_field_destroy(field); + grid_destroy(g); +} + /* ============================================================================ * MAIN * ============================================================================ */ @@ -553,5 +679,6 @@ int main(void) { RUN_TEST(test_heat_source); RUN_TEST(test_thermal_bc_dirichlet_neumann); RUN_TEST(test_thermal_bc_all_periodic_noop); + RUN_TEST(test_thermal_bc_3d_front_back); return UNITY_END(); } From 869b01678366f9ea8fb1e19748efe76839f4422b Mon Sep 17 00:00:00 2001 From: shaia Date: Tue, 7 Apr 2026 20:54:59 +0300 Subject: [PATCH 4/5] Remove unused variable in 3D thermal BC test --- tests/solvers/energy/test_energy_solver.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/solvers/energy/test_energy_solver.c b/tests/solvers/energy/test_energy_solver.c index 494e553..2f4a6b6 100644 --- a/tests/solvers/energy/test_energy_solver.c +++ b/tests/solvers/energy/test_energy_solver.c @@ -563,8 +563,6 @@ static void test_thermal_bc_3d_front_back(void) { flow_field* field = flow_field_create(TBC3D_NX, TBC3D_NY, TBC3D_NZ); TEST_ASSERT_NOT_NULL(field); - size_t total = TBC3D_PLANE * TBC3D_NZ; - /* Fill T with a recognizable 3D pattern: T = k*1000 + j*100 + i */ for (size_t k = 0; k < TBC3D_NZ; k++) { for (size_t j = 0; j < TBC3D_NY; j++) { From a9dd20c1d37fa5a2d6e2124eb85dbb529819a83c Mon Sep 17 00:00:00 2001 From: shaia Date: Tue, 7 Apr 2026 20:57:21 +0300 Subject: [PATCH 5/5] Fix test assertions for periodic and 3D BC application order Update all-periodic test to verify periodic copy behavior instead of no-op (periodic now actively copies from opposite interior cells). Fix 3D test: back/front faces run after bottom/top, so z-face BCs overwrite corner intersections at k=0 and k=nz-1. --- tests/solvers/energy/test_energy_solver.c | 98 ++++++++++++++++------- 1 file changed, 68 insertions(+), 30 deletions(-) diff --git a/tests/solvers/energy/test_energy_solver.c b/tests/solvers/energy/test_energy_solver.c index 2f4a6b6..92c83c2 100644 --- a/tests/solvers/energy/test_energy_solver.c +++ b/tests/solvers/energy/test_energy_solver.c @@ -502,37 +502,80 @@ static void test_thermal_bc_dirichlet_neumann(void) { } /* ============================================================================ - * TEST 8: Thermal BC — All Periodic (no-op) + * TEST 8: Thermal BC — All Periodic * - * Default thermal_bc config (all PERIODIC) should not modify T at all. + * Default thermal_bc config (all PERIODIC) copies boundary cells from the + * opposite interior cell. Verify: + * - Interior cells are unchanged + * - Boundary cells equal the opposite interior cell * ============================================================================ */ -static void test_thermal_bc_all_periodic_noop(void) { - grid* g = grid_create(9, 9, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); +#define PER_NX 9 +#define PER_NY 9 + +static void test_thermal_bc_all_periodic(void) { + grid* g = grid_create(PER_NX, PER_NY, 1, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); TEST_ASSERT_NOT_NULL(g); grid_initialize_uniform(g); - flow_field* field = flow_field_create(9, 9, 1); + flow_field* field = flow_field_create(PER_NX, PER_NY, 1); TEST_ASSERT_NOT_NULL(field); - /* Fill T with distinct values */ - for (size_t n = 0; n < 81; n++) { - field->T[n] = 1000.0 + (double)n; - field->rho[n] = 1.0; + /* Fill T with a recognizable pattern: T = 100 + i + j*10 */ + for (size_t j = 0; j < PER_NY; j++) { + for (size_t i = 0; i < PER_NX; i++) { + size_t idx = j * PER_NX + i; + field->T[idx] = 100.0 + (double)i + (double)j * 10.0; + field->rho[idx] = 1.0; + } } - /* Save original T */ - double T_orig[81]; - memcpy(T_orig, field->T, 81 * sizeof(double)); + /* Save original interior T for comparison */ + double T_orig[PER_NX * PER_NY]; + memcpy(T_orig, field->T, sizeof(T_orig)); ns_solver_params_t params = ns_solver_params_default(); - params.alpha = 0.01; /* Energy enabled, but thermal_bc is all-periodic */ + params.alpha = 0.01; /* Energy enabled, all-periodic by default */ energy_apply_thermal_bcs(field, ¶ms); - /* T should be completely unchanged */ - for (size_t n = 0; n < 81; n++) { - TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[n], field->T[n]); + /* Interior cells (1..nx-2, 1..ny-2) should be unchanged */ + for (size_t j = 1; j < PER_NY - 1; j++) { + for (size_t i = 1; i < PER_NX - 1; i++) { + size_t idx = j * PER_NX + i; + TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[idx], field->T[idx]); + } + } + + /* Left face (i=0): should equal interior column i=nx-2 (before bottom/top overwrite) */ + /* Bottom/top periodic run after left/right, so j=0 and j=ny-1 are overwritten. + * Check interior rows only. */ + for (size_t j = 1; j < PER_NY - 1; j++) { + size_t left_idx = j * PER_NX; + size_t src_idx = j * PER_NX + (PER_NX - 2); + TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[src_idx], field->T[left_idx]); + } + + /* Right face (i=nx-1): should equal interior column i=1 */ + for (size_t j = 1; j < PER_NY - 1; j++) { + size_t right_idx = j * PER_NX + (PER_NX - 1); + size_t src_idx = j * PER_NX + 1; + TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[src_idx], field->T[right_idx]); + } + + /* Bottom face (j=0): copies from j=ny-2, but left/right periodic may have + * modified i=0 and i=nx-1 in that row. Check interior columns. */ + for (size_t i = 1; i < PER_NX - 1; i++) { + size_t bot_idx = i; + size_t src_idx = (PER_NY - 2) * PER_NX + i; + TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[src_idx], field->T[bot_idx]); + } + + /* Top face (j=ny-1): copies from j=1 */ + for (size_t i = 1; i < PER_NX - 1; i++) { + size_t top_idx = (PER_NY - 1) * PER_NX + i; + size_t src_idx = PER_NX + i; + TEST_ASSERT_DOUBLE_WITHIN(1e-15, T_orig[src_idx], field->T[top_idx]); } flow_field_destroy(field); @@ -591,18 +634,12 @@ static void test_thermal_bc_3d_front_back(void) { energy_apply_thermal_bcs(field, ¶ms); - /* Verify back face (k=0): all cells should be 50.0 (Dirichlet) */ + /* Verify back face (k=0): back Dirichlet runs after bottom/top, + * so the entire k=0 plane should be 50.0. */ for (size_t j = 0; j < TBC3D_NY; j++) { for (size_t i = 0; i < TBC3D_NX; i++) { size_t idx = j * TBC3D_NX + i; - /* Bottom/top Dirichlet overwrites corners on j=0 and j=ny-1 */ - if (j == 0) { - TEST_ASSERT_DOUBLE_WITHIN(1e-15, 99.0, field->T[idx]); - } else if (j == TBC3D_NY - 1) { - TEST_ASSERT_DOUBLE_WITHIN(1e-15, 77.0, field->T[idx]); - } else { - TEST_ASSERT_DOUBLE_WITHIN(1e-15, 50.0, field->T[idx]); - } + TEST_ASSERT_DOUBLE_WITHIN(1e-15, 50.0, field->T[idx]); } } @@ -642,16 +679,17 @@ static void test_thermal_bc_3d_front_back(void) { } } - /* Verify bottom face (j=0): Dirichlet 99.0 across all k */ - for (size_t k = 0; k < TBC3D_NZ; k++) { + /* Verify bottom face (j=0): Dirichlet 99.0 at interior k planes. + * k=0 is overwritten by back Dirichlet, k=nz-1 by front Neumann. */ + for (size_t k = 1; k < TBC3D_NZ - 1; k++) { size_t base = k * TBC3D_PLANE; for (size_t i = 0; i < TBC3D_NX; i++) { TEST_ASSERT_DOUBLE_WITHIN(1e-15, 99.0, field->T[base + i]); } } - /* Verify top face (j=ny-1): Dirichlet 77.0 across all k */ - for (size_t k = 0; k < TBC3D_NZ; k++) { + /* Verify top face (j=ny-1): Dirichlet 77.0 at interior k planes */ + for (size_t k = 1; k < TBC3D_NZ - 1; k++) { size_t base = k * TBC3D_PLANE; for (size_t i = 0; i < TBC3D_NX; i++) { TEST_ASSERT_DOUBLE_WITHIN(1e-15, 77.0, @@ -676,7 +714,7 @@ int main(void) { RUN_TEST(test_backward_compatibility); RUN_TEST(test_heat_source); RUN_TEST(test_thermal_bc_dirichlet_neumann); - RUN_TEST(test_thermal_bc_all_periodic_noop); + RUN_TEST(test_thermal_bc_all_periodic); RUN_TEST(test_thermal_bc_3d_front_back); return UNITY_END(); }