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..663d299 100644 --- a/lib/include/cfd/solvers/energy_solver.h +++ b/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 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..5303746 100644 --- a/lib/src/solvers/energy/cpu/energy_solver.c +++ b/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 = ¶ms->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]; + } + } + } +} 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..92c83c2 100644 --- a/tests/solvers/energy/test_energy_solver.c +++ b/tests/solvers/energy/test_energy_solver.c @@ -439,6 +439,268 @@ 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 + * + * 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 + * ============================================================================ */ + +#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(PER_NX, PER_NY, 1); + TEST_ASSERT_NOT_NULL(field); + + /* 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 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, all-periodic by default */ + + energy_apply_thermal_bcs(field, ¶ms); + + /* 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); + 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); + + /* 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): 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; + 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 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 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, + field->T[base + (TBC3D_NY - 1) * TBC3D_NX + i]); + } + } + + flow_field_destroy(field); + grid_destroy(g); +} + /* ============================================================================ * MAIN * ============================================================================ */ @@ -451,5 +713,8 @@ 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); + RUN_TEST(test_thermal_bc_3d_front_back); 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 ---- */