SPH Solver Using Gaussian Smoothing Kernel

sphluids.cpp

#define BOOST_LOG_DYN_LINK 1

#include <iostream>
#include <boost/log/trivial.hpp>
#include <SDL2/SDL.h>
#include <sys/time.h>
#include <cmath>

using namespace std;

const int SCREEN_WIDTH = 1024;
const int SCREEN_HEIGHT = 512;
const int PARTICLE_WIDTH = 10;
const int PARTICLE_COUNT = 60;
const int SDL_DELAY = 8;
const double CONST_H = 0.01;
const double CONST_K = 0.1;
const double CONST_N = 1;
const double MAX_FORCE = 0.8;
const int METABALL_WIDTH = 5;

const double particleMass = 1;
bool mouseDown = false;

struct Vec3 {
    double x, y, z;

    Vec3(double a=0, double b=0, double c=0) : x(a), y(b), z(c) {}

    void add(double a, double b, double c) {
        x += a; y += b; z += c;
    }

    inline double dist() const {
        return sqrt(x*x+y*y+z*z);
    }

    Vec3& operator+=(const Vec3& other) {
        x += other.x; y += other.y, z+= other.z;
        return *this;
    }

    Vec3 operator+(const Vec3& other) const {
        return Vec3(x+other.x, y+other.y, z+other.z);
    }

    Vec3 operator-(const Vec3& other) const {
        return Vec3(x-other.x, y-other.y, z-other.z);
    }
};

Vec3 operator*(const Vec3& v, double other) { return Vec3(v.x*other, v.y*other, v.z*other); }
Vec3 operator*(double other, const Vec3& v) { return Vec3(v.x*other, v.y*other, v.z*other); }
Vec3 operator/(const Vec3& v, double other) { 
    other = 1/other;
    return Vec3(v.x*other, v.y*other, v.z*other); 
}

ostream& operator<<(ostream& os, const Vec3& v) {
    os << '(' << v.x << ", " << v.y << ", " << v.z << ')';
    return os;
}

const Vec3 g{0, -0.1, 0};

struct App {
    SDL_Window *window;
    SDL_Renderer *renderer;
};

struct Particle {
    SDL_Rect rect;

    Vec3 force;
    Vec3 velocity;
    Vec3 acc_prev;

    // position is normalized [0, 1]
    Vec3 pos;

    double mass = particleMass;

    Particle() {
        rect.w = PARTICLE_WIDTH; 
        rect.h = PARTICLE_WIDTH;
    }

    Particle(double x, double y) {
        pos.x = x; pos.y = y; pos.z = 0.0;
        rect.x = x*SCREEN_WIDTH-PARTICLE_WIDTH/2;
        rect.y = SCREEN_HEIGHT-y*SCREEN_HEIGHT-PARTICLE_WIDTH/2;
        rect.w = PARTICLE_WIDTH; rect.h = PARTICLE_WIDTH;
    }

    inline void setPos(double x, double y, double z) {
        pos.x = x; pos.y = y; pos.z = z;
        rect.x = x*SCREEN_WIDTH-PARTICLE_WIDTH/2;
        rect.y = SCREEN_HEIGHT-y*SCREEN_HEIGHT-PARTICLE_WIDTH/2;
    }

    inline void setPos(const Vec3& v) {
        pos.x = v.x; pos.y = v.y; pos.z = v.z;
        rect.x = v.x*SCREEN_WIDTH-PARTICLE_WIDTH/2;
        rect.y = SCREEN_HEIGHT-v.y*SCREEN_HEIGHT-PARTICLE_WIDTH/2;
    }

    void tick(double dt) {
        /* https://en.wikipedia.org/wiki/Leapfrog_integration#Algorithm */
        Vec3 acc = force/mass;
        pos = pos + velocity*dt + 0.5*acc*dt*dt;
        velocity = velocity + 0.5*(acc_prev + acc)*dt;
        acc_prev = acc;

        setPos(pos);
    }
};

App app;
Particle particles[PARTICLE_COUNT];

union SDL_Event;

/* Why is there no straightforward time() https://stackoverflow.com/a/44896326 */
long long time() {
    struct timeval tv;

    gettimeofday(&tv,NULL);
    return (((long long)tv.tv_sec)*1000)+(tv.tv_usec/1000);
}

/*  Weighting/Smoothing function (Gaussian Kernel)
    https://github.com/pmocz/sph-python/blob/master/sph.py 
*/
double W(const Vec3& vec) {
    // h is smoothing length.
    static const double h2_inv = 1/(CONST_H*CONST_H);
    static const double w_coef = pow(1.0/(CONST_H*sqrt(M_PI)), 3.0);

    double r = vec.dist();
    double w = w_coef * exp(-r*r*h2_inv);

    // pairwise w_ij
    return w;
}

/* https://github.com/pmocz/sph-python/blob/master/sph.py */
Vec3 gradW(const Vec3& vec) {
    // h is smoothing length.
    static const double h2 = CONST_H*CONST_H;
    static const double h5 = h2*h2*CONST_H;
    static const double pi_15 = pow(M_PI, 1.5);
    static const double denom_inv = 1/(h5*pi_15);

    double r = vec.dist();
    double n = -2.0*exp(-r*r/(h2))*denom_inv;

    // w_x, w_y, w_z
    return Vec3(n*vec.x, n*vec.y, n*vec.z);
}

double density(const Particle& p, int len) {
    double rho = 0.0;

    for (int i = 0; i < len; i++) {
        rho += particles[i].mass*W(particles[i].pos-p.pos);
    }

    return rho;
}

double pressure(double rho, double k, double n) {
    return k * pow(rho, 1+1/n);
}

Vec3 fluidAccel(const Particle& p, double h, double k, double n) {
    double rho = density(p, PARTICLE_COUNT);
    double P = pressure(rho, k, n);

    double ax, ay, az; ax = ay = az = 0.0;

    for (int i = 0; i < PARTICLE_COUNT; i++) {
        Vec3 dpos = particles[i].pos - p.pos;
        Vec3 dW = gradW(dpos);

        double rho_j = density(particles[i], PARTICLE_COUNT);
        double P_j = pressure(rho_j, k, n);
        double coef = 0.000005*particles[i].mass*(P/(rho*rho)+P_j/(rho_j*rho_j));
        ax += coef*dW.x;
        ay += coef*dW.y;
        az += coef*dW.z;
    }

    return Vec3(ax, ay, az) - 0.8*p.velocity;
}

void doInput() {
    SDL_Event event;

    while (SDL_PollEvent(&event)) {
        switch (event.type) {
            case SDL_QUIT:
                exit(0);
                break;
            case SDL_MOUSEBUTTONDOWN:
                mouseDown = true;
                break;
            case SDL_MOUSEBUTTONUP:
                mouseDown = false;
                break;
            default:
                break;
        }
    }
}

void initSimulation() {
    for (int i = 0; i < PARTICLE_COUNT; i++) {
        particles[i].setPos(rand()/(double)RAND_MAX, rand()/(double)RAND_MAX, 0.0);
        particles[i].velocity = {rand()/(double)RAND_MAX-0.5, rand()/(double)RAND_MAX-0.5, 0.0};
    }
}

void doPhysics(double dt) {
    int x, y;
    SDL_GetMouseState(&x, &y);
    Vec3 cur(x/(double)SCREEN_WIDTH, 1-y/(double) SCREEN_HEIGHT);

    for (int i = 0; i < PARTICLE_COUNT; i++) {
        if (particles[i].pos.y < 0.0 || particles[i].pos.y > 1) {
            if (particles[i].pos.y < 0) particles[i].pos.y = 0.01;
            else particles[i].pos.y = 0.95;
            particles[i].velocity.y *= -0.95;
        }
        if (particles[i].pos.x < 0.0 || particles[i].pos.x > 1) {
            if (particles[i].pos.x < 0) particles[i].pos.x = 0.01;
            else particles[i].pos.x = 0.99;
            particles[i].velocity.x *= -0.95;
        }

        particles[i].force = fluidAccel(particles[i], CONST_H, CONST_K, CONST_N);
        
        double mag_F = particles[i].force.dist();
        if (mag_F > MAX_FORCE) 
            particles[i].force = particles[i].force/mag_F * MAX_FORCE;
        
        if (mouseDown) {
            particles[i].force += cur-particles[i].pos;
        } else {
            particles[i].force += g;
        }

        particles[i].tick(dt);
    }
}

void prepareScene() {
    SDL_SetRenderDrawColor(app.renderer, 100, 100, 100, 255);
    SDL_RenderClear(app.renderer);

    SDL_SetRenderDrawColor(app.renderer, 150, 150, 250, 255);
    SDL_Rect blob_rect{0, 0, METABALL_WIDTH, METABALL_WIDTH};

    static const int ball_offset = METABALL_WIDTH/2;

    for (int i = 0; i < SCREEN_WIDTH; i+=METABALL_WIDTH) {
        double x = i/(double)SCREEN_WIDTH;
        blob_rect.x = i-ball_offset;
        for (int j = 0; j < SCREEN_HEIGHT; j+=METABALL_WIDTH) {
            blob_rect.y = j+ball_offset;
            double y = 1-j/(double)SCREEN_HEIGHT;
            double metaball_v = 0.0;
            for (int i = 0; i < PARTICLE_COUNT; i++) {
                double r_x = (particles[i].pos.x-x);
                double r_y = (particles[i].pos.y-y);
                metaball_v += 1/(r_x*r_x + r_y*r_y);
            }

            if (metaball_v > 1500) {
                SDL_RenderFillRect(app.renderer, &blob_rect);
            }
        }
    }

    SDL_SetRenderDrawColor(app.renderer, 170, 170, 255, 255);
    for (int i = 0; i < PARTICLE_COUNT; i++) {
        SDL_RenderFillRect(app.renderer, &particles[i].rect);
    }
}

void presentScene() {
    SDL_RenderPresent(app.renderer);
}

void initSDL() {
    BOOST_LOG_TRIVIAL(info) << "Starting SDL...\n";

    if(SDL_Init(SDL_INIT_VIDEO|SDL_INIT_AUDIO)==-1) { 
        BOOST_LOG_TRIVIAL(fatal) << "Couldn't initialize video or audio\n";
        exit(-1);
    }

    BOOST_LOG_TRIVIAL(info) << "Initialized SDL\n";

    app.window = SDL_CreateWindow("SPHluids", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED,
        SCREEN_WIDTH, SCREEN_HEIGHT, 0);

    if (!app.window) {
        BOOST_LOG_TRIVIAL(fatal) << "Couldn't initialize window\n";
    }

    SDL_SetHint(SDL_HINT_RENDER_SCALE_QUALITY, "linear");

    app.renderer = SDL_CreateRenderer(app.window, -1, SDL_RENDERER_ACCELERATED);

    if (!app.renderer) {
        BOOST_LOG_TRIVIAL(fatal) << "Couldn't initialize renderer\n";
    }
}

int main() {
    initSDL();
    initSimulation();

    double old = time();
    while (true) {
        double now = time();
        double dt = now-old;
        doPhysics(dt*0.005);
        old = now;

        prepareScene();
        doInput();
        presentScene();
        
        SDL_SetWindowTitle(app.window, to_string(1000/dt).c_str());
        SDL_Delay(SDL_DELAY);
    }

    return 0;
}