a

CMakeLists.txt

@@ -9,11 +9,24 @@ find_package(SDL2 REQUIRED)
 
 add_subdirectory(glm)
 
-set(CMAKE_CUDA_FLAGS ${CMAKE_CUDA_FLAGS} "-G")
+set(CMAKE_CUDA_FLAGS ${CMAKE_CUDA_FLAGS} "")
 
-add_executable(cudagravity main.cu particle.cu particle.cuh definitions.cuh cudaMemory.cuh sdlHelper.cu sdlHelper.cuh)
+# Simulator
 
-target_link_libraries(cudagravity PUBLIC glm::glm SDL2::SDL2)
+add_executable(cudagravity_sim sim/main.cu common/particle.hpp common/definitions.hpp sim/cudaMemory.cuh sim/exportHelper.cu sim/exportHelper.cuh)
 
-set_target_properties(cudagravity PROPERTIES
+target_include_directories(cudagravity_sim PUBLIC common)
+
+target_link_libraries(cudagravity_sim PUBLIC glm::glm)
+
+set_target_properties(cudagravity_sim PROPERTIES
         CUDA_SEPARABLE_COMPILATION ON)
+
+# Viewer
+
+find_package(PkgConfig)
+pkg_check_modules(SDL2GFX REQUIRED SDL2_gfx)
+
+add_executable(cudagravity_viewer viewer/main.cpp common/particle.hpp common/definitions.hpp)
+target_include_directories(cudagravity_viewer PUBLIC common ${SDL2GFX_INCLUDE_DIRS}})
+target_link_libraries(cudagravity_viewer PUBLIC glm::glm SDL2::SDL2 ${SDL2GFX_LIBRARIES})

common/definitions.hpp

@@ -0,0 +1,3 @@
+#include "glm/vec2.hpp"
+
+typedef glm::vec<2, float, glm::defaultp> vec2;

common/particle.hpp

@@ -0,0 +1,14 @@
+#ifndef VERLET_PARTICLE_HPP
+#define VERLET_PARTICLE_HPP
+
+#include "definitions.hpp"
+
+class Particle {
+public:
+    float mass;
+    vec2 position;
+    vec2 velocity;
+};
+
+
+#endif //VERLET_PARTICLE_HPP

definitions.cuh

@@ -1,4 +0,0 @@
-#include <glm/vec2.hpp>
-
-typedef glm::vec<2, float, glm::defaultp> vec2;
-#define CUDA_CALL(CALL) if (CALL != cudaSuccess) { printf("CUDA call failed at %s:%i\n", __FILE__,__LINE__); exit(-1); }

particle.cu

@@ -1,5 +0,0 @@
-//
-// Created by kai on 02/02/23.
-//
-
-#include "particle.cuh"

particle.cuh

@@ -1,16 +0,0 @@
-
-
-#ifndef VERLET_PARTICLE_CUH
-#define VERLET_PARTICLE_CUH
-
-#include "definitions.cuh"
-
-class Particle {
-public:
-    float mass;
-    vec2 position;
-    vec2 velocity;
-};
-
-
-#endif //VERLET_PARTICLE_CUH

sdlHelper.cu

@@ -1,88 +0,0 @@
-//
-// Created by kai on 05/02/23.
-//
-
-#include "sdlHelper.cuh"
-
-sdlHelper::sdlHelper() {
-    if (getenv("NO_DISPLAY") != nullptr) {
-        this->displaying = false;
-        return;
-    }
-
-    if (instance == nullptr) {
-        instance = this;
-
-        this->renderThread = std::thread(&sdlHelper::renderThreadFunc, this);
-    }
-}
-
-sdlHelper::~sdlHelper() {
-    if (instance == this) {
-        running = false;
-
-        this->renderThread.join();
-        SDL_Quit();
-    }
-}
-
-bool sdlHelper::drawParticles(Particle* particles, size_t particleCount) {
-    if (!this->displaying) return this->running;
-    this->particleRects.clear();
-    this->particleRects.reserve(particleCount);
-
-    float maxDistance = 0;
-
-    for (int i = 0; i < particleCount; i++) {
-        maxDistance = std::max(maxDistance, particles[i].position.x);
-        maxDistance = std::max(maxDistance, particles[i].position.y);
-    }
-
-    for (int i = 0; i < particleCount; i++) {
-        SDL_Rect rect;
-
-        vec2 pos = particles[i].position;
-        pos /= (maxDistance * 1.1);
-        pos *= 512;
-        pos += vec2(512, 512);
-
-        rect.x = pos.x - 1;
-        rect.y = pos.y - 1;
-        rect.w = 2;
-        rect.h = 2;
-        this->particleRects.push_back(rect);
-    }
-    return this->running;
-}
-
-sdlHelper *sdlHelper::getInstance() {
-    return this->instance;
-}
-
-void sdlHelper::renderThreadFunc() {
-    this->window = SDL_CreateWindow("cudaGravity", SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED, 1024, 1024, SDL_WINDOW_VULKAN);
-    this->renderer = SDL_CreateRenderer(this->window, -1, SDL_RENDERER_ACCELERATED | SDL_RENDERER_PRESENTVSYNC);
-
-    while (running) {
-        SDL_Event e;
-
-        while (SDL_PollEvent(&e)) {
-            if (e.type == SDL_QUIT) {
-                running = false;
-            }
-        }
-
-        SDL_SetRenderDrawColor(this->renderer, 0, 0, 0, 255);
-        SDL_RenderClear(this->renderer);
-        SDL_SetRenderDrawColor(this->renderer, 255, 255, 255, 255);
-
-        for (auto & particleRect : particleRects) {
-            SDL_RenderFillRect(this->renderer, &particleRect);
-        }
-
-        SDL_RenderPresent(this->renderer);
-    }
-
-    SDL_DestroyRenderer(this->renderer);
-    SDL_DestroyWindow(this->window);
-}

sdlHelper.cuh

@@ -1,35 +0,0 @@
-//
-// Created by kai on 05/02/23.
-//
-
-#ifndef VERLET_SDLHELPER_CUH
-#define VERLET_SDLHELPER_CUH
-
-#include <SDL2/SDL.h>
-#include <vector>
-#include <thread>
-#include "particle.cuh"
-
-class sdlHelper {
-public:
-    sdlHelper();
-    ~sdlHelper();
-
-    bool drawParticles(Particle* particles, size_t particleCount);
-    sdlHelper* getInstance();
-private:
-    void renderThreadFunc();
-
-    sdlHelper* instance = nullptr;
-    SDL_Window* window = nullptr;
-    SDL_Renderer* renderer = nullptr;
-    bool running = true;
-    bool displaying = true;
-
-    std::thread renderThread;
-
-    std::vector<SDL_Rect> particleRects;
-};
-
-
-#endif //VERLET_SDLHELPER_CUH

sim/cudaMemory.cuh

@@ -7,9 +7,10 @@
 #ifndef VERLET_CUDAMEMORY_CUH
 #define VERLET_CUDAMEMORY_CUH
 
-#include "definitions.cuh"
+#include "../common/definitions.hpp"
 #include <cstdio>
 
+#define CUDA_CALL(CALL) if (CALL != cudaSuccess) { printf("CUDA call failed at %s:%i\n", __FILE__,__LINE__); exit(-1); }
 
 template<typename T>
 class CudaMemory {
@@ -46,13 +47,25 @@ public:
     }
 
     void sync() {
-        CUDA_CALL(cudaMemcpy(this->hostPtr, this->devPtr, sizeof(T) * this->N, cudaMemcpyDeviceToHost));
+        CUDA_CALL(cudaMemcpyAsync(this->hostPtr, this->devPtr, sizeof(T) * this->N, cudaMemcpyDeviceToHost));
+    }
+
+    void sync_to_async(T** ptr) {
+        CUDA_CALL(cudaMemcpyAsync(*ptr, this->devPtr, sizeof(T) * this->N, cudaMemcpyDeviceToHost));
+    }
+
+    void sync_to(T** ptr) {
+        CUDA_CALL(cudaMemcpy(*ptr, this->devPtr, sizeof(T) * this->N, cudaMemcpyDeviceToHost));
     }
 
     T& operator[](size_t idx) {
         return this->hostPtr[idx];
     }
 
+    size_t size() {
+        return this->N * sizeof(T);
+    }
+
 private:
     T* hostPtr;
     T* devPtr;

sim/exportHelper.cu

@@ -0,0 +1,104 @@
+//
+// Created by kai on 05/02/23.
+//
+
+#include "exportHelper.cuh"
+
+#include <sys/stat.h>
+
+exportHelper::exportHelper() {
+    mkdir("output", S_IRWXU);
+
+    char* epochs_str = getenv("EPOCHS");
+
+    if (epochs_str != nullptr) {
+        if (sscanf(epochs_str, "%zu", &this->n) == 0) {
+            printf("Failed to parse EPOCHS environment variable - required when NO_DISPLAY is set\n");
+            should_stop = false;
+        }
+    } else {
+        printf("Failed to parse EPOCHS environment variable - required when NO_DISPLAY is set\n");
+        should_stop = false;
+    }
+
+    if (instance == nullptr) {
+        instance = this;
+    }
+}
+
+void exportHelper::stop() {
+    if (instance == this) {
+        should_stop = false;
+
+        for (auto& thread : exportThreads) {
+            if (thread.running) {
+                thread.thread.join();
+            }
+        }
+
+    }
+}
+
+exportHelper::~exportHelper() {
+    stop();
+}
+
+exportHelper *exportHelper::getInstance() {
+    return this->instance;
+}
+
+void exportHelper::imageWriter(size_t idx, thread* thread) {
+    char filename[256];
+    sprintf(filename, "output/%zu.dat", idx);
+
+    FILE* f = fopen64(filename, "wb");
+
+    auto* particles_copy = (Particle*)malloc(this->particles->size());
+    this->particles->sync_to(&particles_copy);
+
+    char data[32 * 5 + 1];
+    data[32*5] = '\0';
+
+    for (size_t i = 0; i < particle_count; i++) {
+	memcpy(data, &particles_copy[i].position.x, sizeof(float));
+        memcpy(&data[32], &particles_copy[i].position.y, sizeof(float));
+        memcpy(&data[64], &particles_copy[i].velocity.x, sizeof(float));
+        memcpy(&data[96], &particles_copy[i].velocity.y, sizeof(float));
+        memcpy(&data[128], &particles_copy[i].mass, sizeof(float));
+
+
+        fprintf(f, "%s\n", particles_copy[i].position.x, particles_copy[i].position.y, particles_copy[i].velocity.x, particles_copy[i].velocity.y, particles_copy[i].mass);
+    }
+
+    free(particles_copy);
+    fflush(f);
+    fclose(f);
+
+    thread->running = false;
+}
+
+void exportHelper::epoch() {
+    this->count++;
+
+    if (this->count % 100 == 0) {
+        printf("Epoch: %zu\n", this->count);
+    }
+
+    if (this->n != 0 && this->count >= this->n) {
+        this->should_stop = false;
+    }
+
+    this->exportThreads.emplace_back();
+
+    this->exportThreads.back().thread = std::thread(&exportHelper::imageWriter, this, this->count, &this->exportThreads.back());
+
+    for (auto it = this->exportThreads.begin(); it != this->exportThreads.end();) {
+        if (!it->running) {
+            it->thread.join();
+            it = this->exportThreads.erase(it);
+        } else {
+            it++;
+        }
+    }
+
+}

sim/exportHelper.cuh

@@ -0,0 +1,53 @@
+//
+// Created by kai on 05/02/23.
+//
+
+#ifndef VERLET_EXPORTHELPER_CUH
+#define VERLET_EXPORTHELPER_CUH
+
+#include <vector>
+#include <thread>
+#include <list>
+#include <atomic>
+#include <particle.hpp>
+#include "cudaMemory.cuh"
+
+struct thread {
+    std::atomic<bool> running = true;
+    std::thread thread = std::thread();
+};
+
+class exportHelper {
+public:
+    bool should_stop = true;
+
+    exportHelper();
+    ~exportHelper();
+
+    exportHelper* getInstance();
+
+    void setParticles(CudaMemory<Particle>* _particles, size_t N) {
+        this->particles = _particles;
+        this->particle_count = N;
+    }
+
+    void epoch();
+
+    void stop();
+private:
+    exportHelper* instance = nullptr;
+    bool displaying = true;
+
+    size_t n = 0;
+    size_t count = 0;
+
+    CudaMemory<Particle>* particles;
+    size_t particle_count;
+
+    std::list<thread> exportThreads;
+
+    void imageWriter(size_t idx, thread* thread);
+};
+
+
+#endif //VERLET_EXPORTHELPER_CUH

sim/main.cu

@@ -1,13 +1,13 @@
 #include <iostream>
-#include "definitions.cuh"
-#include "particle.cuh"
+#include "../common/definitions.hpp"
+#include "../common/particle.hpp"
 #include "cudaMemory.cuh"
-#include "sdlHelper.cuh"
-#include <glm/glm.hpp>
-#include <random>
+#include "exportHelper.cuh"
+#include "glm/glm.hpp"
+#include <unistd.h>
 
 __constant__ const float G_CONSTANT = 1;
-__constant__ float dt = 0.1;
+__constant__ float dt = 0.0001;
 __constant__ const float softening = 0.0000;
 __constant__ const dim3 threadsPerBlock_gpu = dim3(1024);
 const dim3 threadsPerBlock_cpu = dim3(1024);
@@ -23,19 +23,14 @@ __global__ void getAccels(Particle* particles, vec2* accels, size_t N, size_t i)
             return;
         }
 
-        Particle particle = particles[i];
-        Particle other_particle = particles[j];
+        Particle p2 = particles[i];
+        Particle p1 = particles[j];
 
-        vec2 d = {
-                other_particle.position.x - particle.position.x,
-                other_particle.position.y - particle.position.y
-        };
+        float dist = glm::distance(p1.position, p2.position);
 
-        float dist = glm::distance(particle.position, other_particle.position);
+        vec2 f = -G_CONSTANT * ((p2.mass * p1.mass) / (dist * dist)) * glm::normalize(p2.position - p1.position);
 
-        vec2 f = -G_CONSTANT * ((particle.mass * other_particle.mass) / (dist * dist)) * glm::normalize(other_particle.position - particle.position);
-
-        vec2 accel = f / particle.mass;
+        vec2 accel = f / p2.mass;
 
         atomicAdd(&(accels[i].x), accel.x);
         atomicAdd(&(accels[i].y), accel.y);
@@ -63,30 +58,25 @@ __global__ void run_step(Particle* particles, vec2* accels, size_t N) {
     }
 }
 
-int main() {
-    sdlHelper sdlHelper;
+float get_random() {
+    float r = arc4random();
 
-    size_t N = 20000;
+    return (r / ((float)UINT32_MAX / 2.0f)) - 1;
+}
+
+int main() {
+    exportHelper exportHelper;
+
+    size_t N = 2000;
 
     CudaMemory<Particle> particles = CudaMemory<Particle>(N);
     CudaMemory<vec2> accelerations = CudaMemory<vec2>(N);
 
-    std::default_random_engine generator(std::random_device{}());
-    std::uniform_real_distribution<float> mass_distribution(1.0, 100000.0);
-    std::uniform_real_distribution<float> pos_distribution(-0.1, 0.1);
-    std::uniform_real_distribution<float> vel_distribution(-2.0, 2.0);
-
-    for (size_t i = 0; i < N; i++) {
-        accelerations[i] = {0, 0};
-        particles[i] = {
-                mass_distribution(generator),
-                {pos_distribution(generator), pos_distribution(generator)},
-                {vel_distribution(generator), vel_distribution(generator)}
-        };
+    for (int i = 0; i < N; i++) {
+        accelerations[i] = vec2(0, 0);
+        particles[i] = Particle{10, vec2(get_random(), get_random()), vec2(get_random(), get_random())};
     }
 
-    cudaDeviceGetAttribute((int*)(&threadsPerBlock_cpu.x), cudaDevAttrMaxThreadsPerBlock, 0);
-
     dim3 n_blocksPerGrid = dim3(
             ((N) + threadsPerBlock_cpu.x - 1) / threadsPerBlock_cpu.x
     );
@@ -100,18 +90,24 @@ int main() {
     particles.send();
     accelerations.send();
 
-    for (int i = 0; i < 2; i++) {
+    exportHelper.setParticles(&particles, N);
+
+    while(true) {
         run_step<<<n_blocksPerGrid, threadsPerBlock_cpu>>>(particles.getDevicePointer(), accelerations.getDevicePointer(), N);
 
         cudaDeviceSynchronize();
-        particles.sync();
+
         accelerations.send();
 
-        if (!sdlHelper.drawParticles(particles.getPointer(), N)) {
+        exportHelper.epoch();
+
+        if (!exportHelper.should_stop) {
             break;
         }
     }
 
+    exportHelper.stop();
+
     cudaDeviceSynchronize();
 
     return 0;

viewer/main.cpp

@@ -0,0 +1,217 @@
+#include <cstdlib>
+#include <SDL2/SDL.h>
+#include <SDL2_rotozoom.h>
+#include <SDL2_framerate.h>
+#include <vector>
+#include <unistd.h>
+#include "definitions.hpp"
+#include "particle.hpp"
+
+
+int main() {
+    double dt = 0.01;
+    double t = 1.0;
+    size_t current_idx = 1;
+    bool pause = true;
+
+    // Toggles
+
+    bool show_fps = false;
+    bool show_help = true;
+    bool show_zoom = true;
+    bool show_state = true;
+    bool show_speed = true;
+    bool show_grid = true;
+    bool auto_size = true;
+    bool shift = false;
+    bool ctrl = false;
+
+    bool update_required = true;
+
+    // Init SDL
+
+    SDL_Init(SDL_INIT_VIDEO);
+
+    SDL_Window* window = SDL_CreateWindow("Verlet", SDL_WINDOWPOS_UNDEFINED, SDL_WINDOWPOS_UNDEFINED, 800, 800, SDL_WINDOW_SHOWN);
+
+    SDL_Event event;
+
+    SDL_Surface* surf = SDL_CreateRGBSurface(0, 800, 800, 32, 0, 0, 0, 0);
+
+    std::vector<Particle> particles;
+
+    FPSmanager fps;
+
+    SDL_initFramerate(&fps);
+    SDL_setFramerate(&fps, 60);
+
+    float max_dim = 0.0;
+
+    while (true) {
+        while (SDL_PollEvent(&event)) {
+            if (event.type == SDL_QUIT) {
+                exit(0);
+            } else if (event.type == SDL_KEYDOWN) {
+                switch (event.key.keysym.sym) {
+                    case SDLK_RSHIFT:
+                    case SDLK_LSHIFT:
+                        shift = true;
+                        break;
+
+                    case SDLK_RCTRL:
+                    case SDLK_LCTRL:
+                        ctrl = true;
+                        break;
+
+                    case SDLK_SPACE:
+                        pause = !pause;
+                        printf("Paused: %s\n", pause ? "true" : "false");
+                        break;
+
+                    case SDLK_LEFT:
+                        t -= dt;
+                        current_idx = (int)t;
+                        update_required = true;
+                        break;
+
+                    case SDLK_RIGHT:
+                        t += dt;
+                        current_idx = (int)t;
+                        update_required = true;
+                        break;
+
+                    case SDLK_KP_PLUS:
+                        dt *= 2;
+                        printf("dt: %f\n", dt);
+                        break;
+
+                    case SDLK_KP_MINUS:
+                        dt /= 2;
+                        printf("dt: %f\n", dt);
+                        break;
+
+                    case SDLK_KP_0:
+                        t = 1;
+                        current_idx = 1;
+                        update_required = true;
+                        break;
+
+                    case SDLK_UP:
+                        if (shift) {
+                            max_dim -= 1;
+                        } else if (ctrl) {
+                            max_dim -= 10;
+                        } else {
+                            max_dim -= 0.1;
+                        }
+
+			max_dim = std::abs(max_dim);
+
+                        update_required = true;
+                        break;
+
+                    case SDLK_DOWN:
+                        if (shift) {
+                            max_dim += 1;
+                        } else if (ctrl) {
+                            max_dim += 10;
+                        } else {
+                            max_dim += 0.1;
+                        }
+                        update_required = true;
+                        break;
+
+                    case SDLK_a:
+                        auto_size = !auto_size;
+                        update_required = true;
+                        printf("Auto size: %s\n", auto_size ? "true" : "false");
+                        break;
+                }
+            } else if (event.type == SDL_KEYUP) {
+                switch (event.key.keysym.sym) {
+                    case SDLK_RSHIFT:
+                    case SDLK_LSHIFT:
+                        shift = false;
+                        break;
+
+                    case SDLK_RCTRL:
+                    case SDLK_LCTRL:
+                        ctrl = false;
+                        break;
+                }
+            }
+        }
+
+        if (!pause) {
+            t += dt;
+            current_idx = (int)t;
+            update_required = true;
+        }
+
+        if (update_required) {
+            char filename[256];
+            sprintf(filename, "./output/%zu.dat", current_idx);
+
+            FILE* f = fopen(filename, "rb");
+
+            if (f == nullptr) {
+                printf("File not found: %s\n", filename);
+                printf("End of simulation\n");
+                pause = true;
+
+            } else {
+		SDL_FillRect(surf, nullptr, SDL_MapRGB(surf->format, 0, 0, 0));
+
+		particles.clear();
+
+                char *cur_line = (char *)malloc(64 * sizeof(char));
+                size_t len = 64;
+
+                while (true) {
+                    if (getline(&cur_line, &len, f) <= 0) {
+                        break;
+                    }
+
+                    float x = *(float*)cur_line;
+		    float y = *(float*)cur_line[32];
+                    float vx = *(float*)cur_line[64];
+                    float vy = *(float*)cur_line[96];
+		    float m = *(float*)cur_line[128];
+
+                    particles.push_back(Particle{m, {x, y}, {vx, vy}});
+                }
+
+                fclose(f);
+
+                if (auto_size) {
+                    for (auto &particle: particles) {
+                        max_dim = std::max(max_dim,
+                                           std::max(std::abs(particle.position.x), std::abs(particle.position.y)));
+                    }
+                }
+
+                for (auto &particle: particles) {
+                    particle.position /= (max_dim * 1.1);
+                    particle.position *= 400.0;
+                    particle.position += vec2(400, 400);
+
+                    SDL_Rect rect = {(int) particle.position.x - 1, (int) particle.position.y - 1, 2, 2};
+
+                    SDL_FillRect(surf, &rect, SDL_MapRGB(surf->format, 255, 255, 255));
+                }
+            }
+
+            update_required = false;
+        }
+
+        SDL_framerateDelay(&fps);
+
+        SDL_Surface* window_surface = SDL_GetWindowSurface(window);
+
+        SDL_BlitSurface(surf, NULL, window_surface, NULL);
+
+        SDL_UpdateWindowSurface(window);
+    }
+
+    return 0;
+}