main.cpp

Go to the documentation of this file.

/// @file main.cpp
/// @brief SPH fluid simulation with raylib visualization
///
/// Dam-break scenario: fluid falls onto a hot sphere (90 degC) and a cold sphere
/// (0 degC). Heat diffuses between particles and exchanges with rigid bodies.
///
/// Controls:
///   SPACE       -- pause / resume
///   R           -- reset to initial state
///   LMB drag    -- spawn cool water particles (20 degC)
///   RMB drag    -- spawn hot water particles (85 degC)
///   +/-         -- increase / decrease simulation speed (substeps per frame)
///   ESC         -- quit
 
#include "raylib.h"
#include "fluid.hpp"
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <string>
 
// Screen / domain constants
static constexpr int   SCR_W    = 1000;
static constexpr int   SCR_H    = 700;
static constexpr float DOM_W    = 1.0f;
static constexpr float DOM_H    = 0.7f;
static constexpr float SCALE    = SCR_W / DOM_W; // 1000 px/m
static constexpr float PRAD_PX  = 4.5f;          // visual particle radius
 
// Coordinate helpers
inline Vector2 sim_to_screen(float sx, float sy) {
    return { sx * SCALE, (DOM_H - sy) * SCALE };
}
 
inline void sim_to_screen(float sx, float sy, float& px, float& py) {
    px = sx * SCALE;
    py = (DOM_H - sy) * SCALE;
}
 
inline void screen_to_sim(float px, float py, float& sx, float& sy) {
    sx = px / SCALE;
    sy = DOM_H - py / SCALE;
}
 
// Temperature -> colour (blue=cold -> cyan -> green -> yellow -> red=hot)
static Color temp_color(float T, float T_min, float T_max) {
    float t = (T - T_min) / (T_max - T_min + 1e-4f);
    t = std::clamp(t, 0.0f, 1.0f);
    // HSV: hue sweeps 240 deg (blue) -> 0 deg (red)
    return ColorFromHSV((1.0f - t) * 240.0f, 1.0f, 0.95f);
}
 
// Initial scene setup
static void setup(physics::FluidSolver& solver) {
    solver.clear();
 
    const float h  = solver.params().h;
    const float dx = 0.8f * h; // particle spacing ~= 0.02 m
 
    // Fluid block (dam-break): left column, ambient temperature
    for (float x = 0.03f; x < 0.53f; x += dx) {
        for (float y = 0.03f; y < 0.43f; y += dx) {
            // Tiny deterministic jitter to break symmetry
            float jx = 0.001f * ((static_cast<int>(x * 1000) % 7) - 3);
            float jy = 0.001f * ((static_cast<int>(y * 1000) % 5) - 2);
            solver.add_particle(x + jx, y + jy, 0.0f, 0.0f, 20.0f);
        }
    }
 
    // Hot sphere  -- right side, lower
    physics::RigidBody hot{};
    hot.x = 0.74f; hot.y = 0.22f;
    hot.radius      = 0.10f;
    hot.temperature = 90.0f;
    hot.mass        = 5000.0f;
    hot.fixed       = true;
    hot.is_heat_source = true;
    solver.add_body(hot);
 
    // Cold sphere  -- upper middle
    physics::RigidBody cold{};
    cold.x = 0.52f; cold.y = 0.57f;
    cold.radius      = 0.08f;
    cold.temperature = 0.0f;
    cold.mass        = 5000.0f;
    cold.fixed       = true;
    cold.is_heat_source = true;
    solver.add_body(cold);
}
 
// Draw the temperature legend bar (bottom-right corner)
static void draw_legend(float T_min, float T_max) {
    const int BAR_X = SCR_W - 130;
    const int BAR_Y = 30;
    const int BAR_W = 20;
    const int BAR_H = 200;
    const int steps = 50;
 
    for (int i = 0; i < steps; ++i) {
        float t = static_cast<float>(i) / steps;
        float T = T_min + t * (T_max - T_min);
        Color c = temp_color(T, T_min, T_max);
        int y = BAR_Y + BAR_H - static_cast<int>(t * BAR_H);
        DrawRectangle(BAR_X, y, BAR_W, BAR_H / steps + 1, c);
    }
    DrawRectangleLines(BAR_X, BAR_Y, BAR_W, BAR_H, WHITE);
 
    // Labels
    DrawText(TextFormat("%.0f degC", T_max), BAR_X + BAR_W + 4, BAR_Y,          16, WHITE);
    DrawText(TextFormat("%.0f degC", (T_min + T_max) * 0.5f),
             BAR_X + BAR_W + 4, BAR_Y + BAR_H / 2,   16, WHITE);
    DrawText(TextFormat("%.0f degC", T_min), BAR_X + BAR_W + 4, BAR_Y + BAR_H - 16, 16, WHITE);
    DrawText("Temp", BAR_X - 2, BAR_Y + BAR_H + 4, 14, GRAY);
}
 
// Main
int main() {
    InitWindow(SCR_W, SCR_H, "SPH Fluid + Heat Transfer");
    SetTargetFPS(60);
 
    physics::FluidParams params;
    params.h    = 0.025f;
    params.rho0 = 1000.0f;
    params.gamma= 7;
    params.c0   = 10.0f;
    params.mu   = 8.0f;
    params.mass = params.rho0 * (0.8f * params.h) * (0.8f * params.h);
    params.gx   = 0.0f;
    params.gy   = -9.81f;
    params.dt   = 0.001f;
    params.xmin = 0.0f; params.xmax = DOM_W;
    params.ymin = 0.0f; params.ymax = DOM_H;
    params.restitution = 0.01f;
    params.alpha_T = 0.005f;
    params.h_conv  = 8.0f;
 
    physics::FluidSolver solver(params);
    setup(solver);
 
    bool paused    = false;
    int  substeps  = 6;   // simulation substeps per rendered frame
    bool add_hot   = false;
 
    while (!WindowShouldClose()) {
        if (IsKeyPressed(KEY_SPACE)) paused = !paused;
        if (IsKeyPressed(KEY_R))     setup(solver);
        if (IsKeyPressed(KEY_EQUAL) || IsKeyPressed(KEY_KP_ADD))
            substeps = std::min(substeps + 1, 20);
        if (IsKeyPressed(KEY_MINUS)  || IsKeyPressed(KEY_KP_SUBTRACT))
            substeps = std::max(substeps - 1, 1);
        if (IsKeyPressed(KEY_G)) {
            for (physics::RigidBody& b : solver.bodies())
                b.fixed = false;
        }
 
        // Spawn particles at mouse position
        if (IsMouseButtonDown(MOUSE_BUTTON_LEFT) || IsMouseButtonDown(MOUSE_BUTTON_RIGHT)) {
            add_hot = IsMouseButtonDown(MOUSE_BUTTON_RIGHT);
            Vector2 mp = GetMousePosition();
            float sx, sy;
            screen_to_sim(mp.x, mp.y, sx, sy);
            const float T = add_hot ? 85.0f : 15.0f;
            // Spawn a small cluster
            for (int k = 0; k < 4; ++k) {
                float ox = 0.01f * ((k % 3) - 1);
                float oy = 0.01f * ((k / 3) - 0);
                float nx = sx + ox, ny = sy + oy;
                if (nx > params.xmin && nx < params.xmax &&
                    ny > params.ymin && ny < params.ymax) {
                    solver.add_particle(nx, ny, 0.0f, 0.0f, T);
                }
            }
        }
 
        if (!paused) {
            for (int s = 0; s < substeps; ++s) {
                solver.step();
                // If heat source: maintain rigid body temperatures
                for (physics::RigidBody& b : solver.bodies()) {
                    if (b.is_heat_source) { /* temperature held fixed already */ }
                }
            }
        }
 
        BeginDrawing();
        ClearBackground(Color{18, 18, 28, 255});
 
        // Draw domain boundary
        DrawRectangleLines(0, 0, SCR_W, SCR_H, Color{60, 60, 80, 255});
 
        const float T_min = solver.min_temp();
        const float T_max = solver.max_temp();
 
        // Draw rigid bodies (filled + outline)
        for (const physics::RigidBody& body : solver.bodies()) {
            Vector2 pos = sim_to_screen(body.x, body.y);
            float   r   = body.radius * SCALE;
            Color   fc  = temp_color(body.temperature, T_min, T_max);
            fc.a = 180;
            DrawCircleV(pos, r, fc);
            DrawCircleLines(static_cast<int>(pos.x), static_cast<int>(pos.y),
                            r, Color{220, 220, 220, 230});
 
            // Label
            const char* label = (body.temperature > 50.0f) ? "HOT" : "COLD";
            DrawText(label,
                     static_cast<int>(pos.x) - 12,
                     static_cast<int>(pos.y) - 8, 14,
                     (body.temperature > 50.0f) ? RED : SKYBLUE);
        }
 
        // Draw particles
        for (const physics::Particle& p : solver.particles()) {
            float px, py;
            sim_to_screen(p.x, p.y, px, py);
            Color c = temp_color(p.temperature, T_min, T_max);
            DrawCircleV({px, py}, PRAD_PX, c);
        }
 
        // HUD
        DrawFPS(10, 10);
        DrawText(TextFormat("Particles: %d", (int)solver.particles().size()),
                 10, 36, 20, WHITE);
        DrawText(TextFormat("Substeps/frame: %d  (sim rate: %.2fx)",
                            substeps, substeps * params.dt * 60.0f),
                 10, 60, 18, LIGHTGRAY);
        DrawText(TextFormat("T range: %.1f degC - %.1f degC", T_min, T_max),
                 10, 84, 18, LIGHTGRAY);
 
        // Controls bar
        DrawRectangle(0, SCR_H - 30, SCR_W, 30, Color{0,0,0,180});
        DrawText("SPACE=pause  R=reset  G=drop bodies  LMB=cold water  RMB=hot water  +/-=speed  ESC=quit",
                 8, SCR_H - 22, 14, Color{180, 180, 180, 255});
 
        // Pause overlay
        if (paused) {
            DrawRectangle(0, 0, SCR_W, SCR_H, Color{0, 0, 0, 100});
            DrawText("PAUSED", SCR_W / 2 - 60, SCR_H / 2 - 20, 40, YELLOW);
        }
 
        // Temperature legend
        draw_legend(T_min, T_max);
 
        EndDrawing();
    }
 
    CloseWindow();
    return 0;
}