feat: Add energy equation solver with Boussinesq buoyancy coupling

CMakeLists.txt

@@ -297,6 +297,10 @@ if(BUILD_TESTS)
     add_executable(test_taylor_green_3d tests/validation/test_taylor_green_3d.c)
     add_executable(test_quiescent_3d tests/validation/test_quiescent_3d.c)
 
+    # Energy equation tests
+    add_executable(test_energy_solver tests/solvers/energy/test_energy_solver.c)
+    add_executable(test_natural_convection tests/validation/test_natural_convection.c)
+
     # Mathematical accuracy tests
     add_executable(test_finite_differences tests/math/test_finite_differences.c)
     add_executable(test_finite_differences_3d tests/math/test_finite_differences_3d.c)
@@ -398,6 +402,8 @@ if(BUILD_TESTS)
     target_link_libraries(test_quiescent_3d PRIVATE CFD::Library unity)
     target_link_libraries(test_poiseuille_flow PRIVATE CFD::Library unity)
     target_link_libraries(test_poiseuille_3d PRIVATE CFD::Library unity)
+    target_link_libraries(test_energy_solver PRIVATE CFD::Library unity)
+    target_link_libraries(test_natural_convection PRIVATE CFD::Library unity)
     target_link_libraries(test_finite_differences PRIVATE unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
     target_include_directories(test_finite_differences PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/lib/include)
     target_link_libraries(test_finite_differences_3d PRIVATE unity $<$<NOT:$<PLATFORM_ID:Windows>>:m>)
@@ -519,6 +525,10 @@ if(BUILD_TESTS)
     add_test(NAME PoiseuilleFlowTest COMMAND test_poiseuille_flow)
     add_test(NAME Poiseuille3DTest COMMAND test_poiseuille_3d)
 
+    # Energy equation tests
+    add_test(NAME EnergySolverTest COMMAND test_energy_solver)
+    add_test(NAME NaturalConvectionTest COMMAND test_natural_convection)
+
     # Mathematical accuracy tests
     add_test(NAME FiniteDifferencesTest COMMAND test_finite_differences)
     add_test(NAME FiniteDifferences3DTest COMMAND test_finite_differences_3d)

lib/CMakeLists.txt

@@ -118,6 +118,8 @@ set(CFD_SCALAR_SOURCES
     src/solvers/navier_stokes/cpu/solver_explicit_euler.c
     src/solvers/navier_stokes/cpu/solver_rk2.c
     src/solvers/navier_stokes/cpu/solver_projection.c
+    # Energy equation solver - scalar
+    src/solvers/energy/cpu/energy_solver.c
     # Linear solvers - scalar implementations only
     src/solvers/linear/cpu/linear_solver_jacobi.c
     src/solvers/linear/cpu/linear_solver_sor.c
@@ -530,6 +532,13 @@ else()
     )
 endif()
 
+# Propagate CFD_ENABLE_OPENMP to consumers (tests, examples)
+# The directory-level add_compile_definitions() only affects targets in lib/,
+# but test executables in the root CMakeLists.txt need this for compile-time guards.
+if(OpenMP_C_FOUND)
+    target_compile_definitions(cfd_library INTERFACE CFD_ENABLE_OPENMP)
+endif()
+
 # Symbol Visibility and Export Header (needed for both builds)
 include(GenerateExportHeader)
 # For INTERFACE libraries, we can't generate export header from it, use cfd_core instead
@@ -574,6 +583,7 @@ set(CFD_LIBRARY_HEADERS
     ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/core/cpu_features.h
     ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/solvers/navier_stokes_solver.h
     ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/solvers/poisson_solver.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/solvers/energy_solver.h
     ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/boundary/boundary_conditions.h
     ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/io/vtk_output.h
     ${CMAKE_CURRENT_SOURCE_DIR}/include/cfd/io/csv_output.h

lib/include/cfd/solvers/energy_solver.h

@@ -0,0 +1,68 @@
+/**
+ * @file energy_solver.h
+ * @brief Energy equation solver for temperature advection-diffusion
+ *
+ * Solves the energy equation for incompressible flow:
+ *   dT/dt + u*nabla(T) = alpha * nabla^2(T) + Q
+ *
+ * where alpha = k/(rho*cp) is thermal diffusivity and Q is a heat source.
+ *
+ * Boussinesq buoyancy coupling adds a buoyancy acceleration to the momentum equations:
+ *   a_buoy = -beta * (T - T_ref) * g
+ *
+ * The energy equation is solved as a post-step after velocity advancement,
+ * using the updated velocity field for advection.
+ */
+
+#ifndef CFD_ENERGY_SOLVER_H
+#define CFD_ENERGY_SOLVER_H
+
+#include "cfd/cfd_export.h"
+#include "cfd/core/cfd_status.h"
+#include "cfd/core/grid.h"
+#include "cfd/solvers/navier_stokes_solver.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * Advance the temperature field by one explicit Euler step.
+ *
+ * Computes: T^{n+1} = T^n + dt * (-u*nabla(T) + alpha*nabla^2(T) + Q)
+ *
+ * Only active when params->alpha > 0. Returns immediately otherwise.
+ *
+ * @param field  Flow field (T is updated in-place, u/v/w are read-only)
+ * @param grid   Computational grid
+ * @param params Solver parameters (alpha, heat_source_func, dt)
+ * @param dt     Time step size
+ * @param time   Current physical time (passed to heat_source_func)
+ * @return CFD_SUCCESS, or CFD_ERROR_DIVERGED if NaN detected in T
+ */
+CFD_LIBRARY_EXPORT cfd_status_t energy_step_explicit(flow_field* field, const grid* grid,
+                                                      const ns_solver_params_t* params,
+                                                      double dt, double time);
+
+/**
+ * Compute Boussinesq buoyancy source terms from the temperature field.
+ *
+ * Adds buoyancy force: source_{u,v,w} += -beta * (T[idx] - T_ref) * g_{x,y,z}
+ *
+ * Only active when params->beta != 0. Does nothing otherwise.
+ *
+ * @param T_local   Temperature at the current grid point
+ * @param params    Solver parameters (beta, T_ref, gravity)
+ * @param source_u  Output: added u-momentum buoyancy source
+ * @param source_v  Output: added v-momentum buoyancy source
+ * @param source_w  Output: added w-momentum buoyancy source
+ */
+CFD_LIBRARY_EXPORT void energy_compute_buoyancy(double T_local, const ns_solver_params_t* params,
+                                                 double* source_u, double* source_v,
+                                                 double* source_w);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* CFD_ENERGY_SOLVER_H */

lib/include/cfd/solvers/navier_stokes_solver.h

@@ -78,6 +78,18 @@ typedef void (*ns_source_func_t)(double x, double y, double z, double t,
                                   double* source_u, double* source_v,
                                   double* source_w);
 
+/**
+ * Custom heat source callback function
+ * @param x X coordinate
+ * @param y Y coordinate
+ * @param z Z coordinate (0.0 for 2D)
+ * @param t Physical time
+ * @param context User-provided context pointer
+ * @return Heat source term Q at the given point
+ */
+typedef double (*ns_heat_source_func_t)(double x, double y, double z, double t,
+                                         void* context);
+
 /**
  * Navier-Stokes solver parameters
  */
@@ -99,6 +111,20 @@ typedef struct {
     /* Custom source term callback (NULL = use default) */
     ns_source_func_t source_func;  /**< Custom source term function pointer */
     void* source_context;           /**< User context for source_func */
+
+    /* Energy equation parameters (alpha > 0 enables energy equation).
+     * NOTE: Energy coupling (temperature advection/diffusion and Boussinesq
+     * buoyancy) is currently only implemented in scalar CPU backends
+     * (explicit_euler, projection, rk2). OMP, AVX2, and GPU backends will
+     * return CFD_ERROR_UNSUPPORTED if alpha > 0 or beta != 0. */
+    double alpha;           /**< Thermal diffusivity k/(rho*cp) [m^2/s], 0 = disabled */
+    double beta;            /**< Thermal expansion coefficient [1/K] (Boussinesq) */
+    double T_ref;           /**< Reference temperature [K] for Boussinesq */
+    double gravity[3];      /**< Gravity vector (gx, gy, gz) [m/s^2] */
+
+    /* 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 */
 } ns_solver_params_t;
 
 
@@ -144,6 +170,7 @@ typedef struct {
     double residual;         /**< Final residual norm */
     double max_velocity;     /**< Maximum velocity magnitude */
     double max_pressure;     /**< Maximum pressure */
+    double max_temperature;  /**< Maximum temperature (when energy equation active) */
     double cfl_number;       /**< Actual CFL number used */
     double elapsed_time_ms;  /**< Wall clock time for solve */
     cfd_status_t status;     /**< Status of the solve */
@@ -338,6 +365,7 @@ static inline ns_solver_stats_t ns_solver_stats_default(void) {
     stats.residual = 0.0;
     stats.max_velocity = 0.0;
     stats.max_pressure = 0.0;
+    stats.max_temperature = 0.0;
     stats.cfl_number = 0.0;
     stats.elapsed_time_ms = 0.0;
     stats.status = CFD_SUCCESS;