Initial commit

.gitignore

@@ -0,0 +1,8 @@
+Cpp_Stuff/.cache/clangd/index
+Cpp_Stuff/build
+Cpp_Stuff/compile_commands.json
+Cpp_Stuff/tags
+Python_Stuff/.vscode
+Python_Stuff/__pycache__
+Python_Stuff/new_venv
+Python_Stuff/venv

Cpp_Stuff/CMakeLists.txt

@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 3.22)
+project(particle_sim)
+# set (CMAKE_EXPORT_COMPILE_COMMANDS ON)
+# set(CMAKE_EXPORT_COMPILE_COMMANDS ON CACHE INTERNAL "")
+add_subdirectory(./src)
+#add_subdirectory(./test)
+

Cpp_Stuff/include/custom_kernels.cuh

@@ -0,0 +1,2 @@
+
+void custom_kernel_process_ray_collision(int size, float * norms, float * direction, float * vel, float * pos, int array_stride, float * line_vec);

Cpp_Stuff/include/hal.hpp

@@ -0,0 +1,3 @@
+#include "vector"
+
+std::vector<float> hal_vector_f32_cross_multiply(std::vector<float> & a, std::vector<float> & b);

Cpp_Stuff/include/sim_renderer.h

@@ -0,0 +1,5 @@
+#include <SDL2/SDL.h>
+#include <SDL2/SDL_video.h>
+
+
+int begin_sim(SDL_Renderer *renderer, SDL_Window *win);

Cpp_Stuff/include/textures.hpp

@@ -0,0 +1,5 @@
+#include <SDL2/SDL.h>
+
+
+void init_textures(SDL_Renderer *in_renderer);
+SDL_Texture * get_texture(int color);

Cpp_Stuff/include/ui_particles.hpp

@@ -0,0 +1,18 @@
+#include <SDL2/SDL.h>
+
+enum Color {blue, green, red, yellow};
+
+class ui_particles
+{
+    private:
+        SDL_Surface *bitmapSurface = NULL;
+        SDL_Texture * bitmapTexture = NULL;
+        SDL_Rect Rectangle;
+        SDL_Renderer *renderer;
+
+    public:
+        ui_particles(Color color, SDL_Renderer *in_renderer);
+        ~ui_particles();
+        void set_location(float x, float y);
+        void draw(void);
+};

Cpp_Stuff/include/universe.hpp

@@ -0,0 +1,39 @@
+#include <vector>
+#include "ui_particles.hpp"
+#include <cublas_v2.h>
+#include <curand.h>
+
+struct data3d{
+    float *data;
+};
+
+class universe
+{
+    private:
+        std::vector<ui_particles> my_particles;
+        int num_of_particles;
+        int array_size;
+        cublasHandle_t cublas_handle;
+        float * device_ray_position;
+        float * device_ray_direction;
+        float * device_position;
+        float * device_velocity;
+        float * device_acceleration;
+        float * device_scratch_padA;
+        float * device_scratch_padB;
+        float * host_position;
+        float * host_velocity;
+        float * host_acceleration;
+        curandGenerator_t gen;
+        void process_collisions_with_other_particles(void);
+        void process_collisions_with_rays(void);
+        void process_particles(void);
+     /* data3d */
+    public:
+        universe(int size);
+        ~universe();
+        void get_positions(float ** in_out_buffer, int * size, int * stride);
+        void process(void);
+};
+
+

Cpp_Stuff/include/util.hpp

@@ -0,0 +1,18 @@
+
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <iostream>
+
+#define CHECK_CUDA_ERROR(val) check((val), #val, __FILE__, __LINE__)
+template <typename T>
+void check(T err, const char* const func, const char* const file,
+           const int line)
+{
+    if (int(err) != int(cudaSuccess))
+    {
+        std::cerr << "CUDA Runtime Error at: " << file << ":" << line << ":" << err 
+                  << std::endl;
+        // We don't exit when we encounter CUDA errors in this example.
+        // std::exit(EXIT_FAILURE);
+    }
+}

Cpp_Stuff/src/CMakeLists.txt

@@ -0,0 +1,31 @@
+# set(CMAKE_CXX_OMPILER "nvcc")
+
+
+find_package(SDL2 REQUIRED)
+# find_package(CUDA REQUIRED)
+find_package(CUDAToolkit REQUIRED)
+
+
+file(GLOB project_cuda_src "**/*.cu")
+file(GLOB project_cpp_src "**/*.cpp")
+# set(project_src ${project_cpp_src} ${project_cuda_src})
+include_directories(${CMAKE_PROJECT_NAME} "${PROJECT_SOURCE_DIR}/../include")
+include_directories(${CMAKE_PROJECT_NAME} ${SDL2_INCLUDE_DIRS})
+#include_directories(${CMAKE_PROJECT_NAME} ${CUDA_INCLUDE_DIRS})
+message(STATUS "${SDL2_LIBRARIES}")
+message(STATUS "${CUDA_LIBRARIES}")
+message(STATUS "${CUDAToolkit_ROOT}")
+enable_language(CUDA)
+
+# add_library(cuda_lib STATIC ${project_cuda_src})
+# add_library(particle_sim STATIC ${project_src})
+
+# add_executable(particle_sim_exe "main.cpp" ${project_src})
+add_executable(particle_sim_exe "main.cpp" ${project_cuda_src} ${project_cpp_src})
+set_target_properties(particle_sim_exe PROPERTIES CXX_STANDARD 17)
+target_compile_options(particle_sim_exe PRIVATE -dlink)
+target_link_libraries( particle_sim_exe PRIVATE ${SDL2_LIBRARIES} )
+target_link_libraries( particle_sim_exe PRIVATE ${SDL2_LIBRARIES} )
+# target_link_libraries( particle_sim_exe PRIVATE -lcublas -lcudrt -lnppc -lnpps -lcurand)
+target_link_libraries( particle_sim_exe PRIVATE CUDA::cublas CUDA::nppc CUDA::npps CUDA::curand)
+set_target_properties( particle_sim_exe PROPERTIES CUDA_SEPARABLE_COMPILATION ON)

Cpp_Stuff/src/main.cpp

@@ -0,0 +1,47 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <string.h>
+#include <stdint.h>
+#include <SDL2/SDL.h>
+#include <SDL2/SDL_video.h>
+#include "sim_renderer.h"
+
+int main( int argc, const char* argv[] )
+{
+    SDL_Renderer *renderer = NULL;
+    SDL_Texture *texture = NULL;
+    SDL_Window *win;
+    SDL_Rect srcR, destR;
+    const char name[] = "Game";
+    if (SDL_Init(SDL_INIT_EVERYTHING) != 0)
+    {
+        SDL_Log("Unable to initialize SDL: %s", SDL_GetError());
+        return 1;
+    }
+
+    win = SDL_CreateWindow(name, 0, 0, 800, 800, SDL_WINDOW_RESIZABLE);
+    if (win == NULL)
+    {
+        // In the case that the window could not be made...
+        printf("Could not create window: %s\n", SDL_GetError());
+        return 1;
+    }
+    renderer = SDL_CreateRenderer(win, -1, SDL_RENDERER_ACCELERATED);
+    SDL_SetRenderDrawColor(renderer, 255, 255, 255, 255);
+    SDL_Rect rectangle;
+
+    rectangle.x = 0;
+    rectangle.y = 0;
+    rectangle.w = 50;
+    rectangle.h = 50;
+    SDL_RenderFillRect(renderer, &rectangle);
+    texture = SDL_CreateTexture(renderer, SDL_PIXELFORMAT_RGBA8888, SDL_TEXTUREACCESS_TARGET, 800, 800);
+
+    (void)begin_sim(renderer, win);
+
+    SDL_DestroyTexture(texture);
+    SDL_DestroyRenderer(renderer);
+    SDL_DestroyWindow(win);
+    return 0;
+}

Cpp_Stuff/src/renderer/sim_renderer.cpp

@@ -0,0 +1,132 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <string.h>
+#include <stdint.h>
+#include <SDL2/SDL.h>
+#include <SDL2/SDL_video.h>
+#include <time.h>
+#include <vector>
+#include "universe.hpp"
+#include "util.hpp"
+#include "textures.hpp"
+
+static clock_t start_time, end_time;
+std::vector<int> test;
+
+static void draw_particles(universe * my_universe, std::vector<ui_particles> *my_particles)
+{
+    float * positions;
+    int stride;
+    int size;
+    my_universe->get_positions(&positions, &size, &stride);
+    for(int i = 0; i < size; i++)
+    {
+        // printf("Test\n");
+        (*my_particles)[i].set_location(positions[i], positions[stride+i]);
+        (*my_particles)[i].draw();
+    }
+}
+
+int begin_sim(SDL_Renderer *renderer, SDL_Window *win)
+{
+    SDL_Rect srcR, destR;
+    int num_of_particles = 10000;
+    bool run = true;
+    SDL_Surface *screen;
+    universe my_universe = universe(num_of_particles);
+    init_textures(renderer);
+    std::vector<ui_particles> my_particles;
+    for(int i = 0; i < num_of_particles; i++)
+    {
+        ui_particles temp_particle = ui_particles(Color(i%4), renderer);
+        my_particles.push_back(temp_particle);
+    }
+    // screen = SDL_SetVideoMode(800, 800, 0, 0);
+    while (run)
+    {
+        clock_t start_t, end_t;
+        SDL_Event event;
+        bool redraw = true;
+        while (SDL_PollEvent(&event))
+        {
+            /* handle your event here */
+            if (event.type == SDL_MOUSEBUTTONDOWN)
+            {
+                if (event.button.button == SDL_BUTTON_MIDDLE) // scroll up
+                {
+                    run = false;
+                }
+                else if ((event.button.button == SDL_BUTTON_LEFT)
+                    || (event.button.button == SDL_BUTTON_RIGHT)) // scroll up
+                    {
+                        // bitmapSurface = SDL_LoadBMP("./sdl_sprites/red.bpm");
+                        SDL_SetRenderDrawColor(renderer, 0x00, 0x00, 0x00, 0x00);
+                        SDL_RenderClear(renderer);
+                        my_universe.process();
+                        start_time = clock();
+                        draw_particles(&my_universe, &my_particles);
+                        end_time = clock();
+
+                        SDL_RenderPresent(renderer);
+                        clock_t t = end_time - start_time;
+                        SDL_Log("No. of clicks %ld clicks (%f seconds) to process the incoming click.\n", t, ((float)t) / CLOCKS_PER_SEC);
+                }
+                else
+                {
+                    /* Do nothing. */
+                }
+                
+            }
+            else if (event.type == SDL_QUIT)
+            {
+                run = false;
+            }
+            else if (event.type == SDL_WINDOWEVENT)
+            {
+                switch (event.window.event)
+                {
+                case SDL_WINDOWEVENT_ENTER:
+                case SDL_WINDOWEVENT_SHOWN:
+                case SDL_WINDOWEVENT_MOVED:
+                case SDL_WINDOWEVENT_EXPOSED:
+                case SDL_WINDOWEVENT_RESTORED:
+                    redraw = true;
+                    break;
+                case SDL_WINDOWEVENT_RESIZED:
+                case SDL_WINDOWEVENT_SIZE_CHANGED:
+                case SDL_WINDOWEVENT_MAXIMIZED:
+                    redraw = true;
+                    SDL_GetWindowSize(win, &destR.w, &destR.h);
+                    // ui_resize(renderer, destR.w, destR.h);
+                    break;
+                }
+            }
+            else
+            {
+                /* Do nothing. */
+            }
+        }
+        if (redraw)
+        {
+            SDL_SetRenderDrawColor(renderer, 0x00, 0x00, 0x00, 0x00);
+            SDL_RenderClear(renderer);
+            start_time = clock();
+            for(int i = 0; i < 100; i++)
+            {
+                my_universe.process();
+            }
+            end_time = clock();
+            draw_particles(&my_universe, &my_particles);
+
+            // ui_redraw_board(renderer);
+            SDL_RenderPresent(renderer);
+
+            clock_t t = end_time - start_time;
+            SDL_Log("No. of clicks %ld clicks (%f seconds).\n",
+                   t, ((float)t) / CLOCKS_PER_SEC);
+        }
+        /* do some other stuff here -- draw your app, etc. */
+    }
+    return 0;
+}

Cpp_Stuff/src/renderer/textures.cpp

@@ -0,0 +1,25 @@
+#include "textures.hpp"
+#include <stdlib.h>
+#include <string.h>
+#include <SDL2/SDL.h>
+#include <iostream>
+
+SDL_Surface *bitmapSurface[4]= {};
+SDL_Texture * bitmapTexture[4] = {};
+static std::string color_paths[] = {"blue", "green", "red", "yellow"};
+static std::string sprite_path = "../sdl_sprites/";
+
+void init_textures(SDL_Renderer *in_renderer)
+{
+    for (size_t color = 0; color < 4; color++)
+    {
+        bitmapSurface[color] = SDL_LoadBMP((sprite_path+color_paths[color]+".bmp").c_str());
+        bitmapTexture[color] = SDL_CreateTextureFromSurface(in_renderer, bitmapSurface[color]);
+    }
+    
+}
+
+SDL_Texture * get_texture(int color)
+{
+    return bitmapTexture[color];
+}

Cpp_Stuff/src/renderer/ui_particles.cpp

@@ -0,0 +1,44 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <string.h>
+#include <stdint.h>
+#include <SDL2/SDL.h>
+#include <SDL2/SDL_video.h>
+#include "ui_particles.hpp"
+#include "textures.hpp"
+
+
+static std::string color_paths[] = {"blue", "green", "red", "yellow"};
+static std::string sprite_path = "../sdl_sprites/";
+
+ui_particles::ui_particles(Color color, SDL_Renderer *in_renderer)
+{
+    renderer = in_renderer;
+    bitmapTexture = get_texture(color);
+    Rectangle.x = 0;
+    Rectangle.y = 0;
+    Rectangle.w = 5;
+    Rectangle.h = 5;
+    // std::cout<<sprite_path<<color_paths[color] << "\n";
+}
+
+ui_particles::~ui_particles()
+{
+}
+
+
+
+void ui_particles::set_location(float x, float y)
+{
+    Rectangle.x = int(x-(Rectangle.w/2));
+    Rectangle.y = int(y-(Rectangle.h/2));
+}
+
+// void ui_particles::draw(SDL_Surface * dest_surface)
+void ui_particles::draw(void)
+{
+    SDL_RenderCopy(renderer, bitmapTexture, NULL, &Rectangle);
+    // SDL_BlitSurface(bitmapSurface, NULL, dest_surface, &Rectangle);
+}
+

Cpp_Stuff/src/simulator/universe.cpp

@@ -0,0 +1,218 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <string.h>
+#include <stdint.h>
+#include <SDL2/SDL.h>
+#include <SDL2/SDL_video.h>
+#include <cuda.h>
+#include <cuda_runtime.h>
+#include "universe.hpp"
+#include "npp.h"
+#include "custom_kernels.cuh"
+#include "util.hpp"
+
+#define DELTA_TIME 0.5
+
+
+void universe::process_collisions_with_rays(void)
+{
+    float temp_directions[][3] = 
+    {
+        {1.0,0.0,0.0},
+        {0.0,1.0,0.0},
+        {1.0,0.0,0.0},
+        {0.0,1.0,0.0}
+    };
+
+    float temp_positions[][3] = 
+    {
+        {0.0,0.0,0.0},
+        {0.0,0.0,0.0},
+        {800.0,800.0,800.0},
+        {800.0,800.0,800.0}
+    };
+
+    float alpha = 1.0;
+    float beta = 0.0;
+
+    float * dev_dot_product_result;
+    float * dev_direction_vec;
+    
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(dev_dot_product_result),
+                long(num_of_particles*sizeof(float))));
+
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(dev_direction_vec),
+                long(num_of_particles*3*sizeof(float))));
+
+    for(int i = 0; i < 4; i++)
+    {
+        //Fill scratchpad A with vec1 = line_position - particle_position
+        for (int j = 0; j < 3; j++)
+        {
+            CHECK_CUDA_ERROR(nppsSubCRev_32f(&device_position[array_size * j],
+                        temp_positions[i][j], &device_scratch_padA[num_of_particles * j],
+                        num_of_particles));
+        }
+        alpha = 1.0;
+        beta = 0.0;
+        CHECK_CUDA_ERROR(cublasSgemm(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N,
+                    num_of_particles, 1, 3, &alpha, device_scratch_padA,
+                    num_of_particles, &device_ray_direction[i * 3], 3, &beta,
+                    dev_dot_product_result, num_of_particles));
+
+        //Use dot product to find projections onto the line directions
+        CHECK_CUDA_ERROR(cublasSgemm(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_T,
+                    num_of_particles, 3, 1, &alpha, dev_dot_product_result,
+                    num_of_particles, &device_ray_direction[i * 3], 3, &beta,
+                    device_scratch_padA, num_of_particles));
+
+        //calc norm
+        CHECK_CUDA_ERROR(nppsSqr_32f(device_scratch_padA, dev_direction_vec,
+                    num_of_particles * 3));
+
+        CHECK_CUDA_ERROR(cublasSaxpy(cublas_handle, num_of_particles, &alpha,
+                    &dev_direction_vec[num_of_particles * 1], 1, 
+                    dev_direction_vec, 1));
+
+        CHECK_CUDA_ERROR(cublasSaxpy(cublas_handle, num_of_particles, &alpha,
+                    &dev_direction_vec[num_of_particles * 2], 1,
+                    dev_direction_vec, 1));
+
+        CHECK_CUDA_ERROR(nppsSqrt_32f_I(dev_direction_vec, num_of_particles));
+
+        //Normalize vector
+        CHECK_CUDA_ERROR(nppsDiv_32f_I(dev_direction_vec, device_scratch_padA,
+                    num_of_particles));
+
+        CHECK_CUDA_ERROR(nppsDiv_32f_I(dev_direction_vec, 
+                    &device_scratch_padA[1 * num_of_particles], num_of_particles));
+
+        CHECK_CUDA_ERROR(nppsDiv_32f_I(dev_direction_vec, 
+                    &device_scratch_padA[2 * num_of_particles], num_of_particles));
+
+        custom_kernel_process_ray_collision(num_of_particles, dev_direction_vec, 
+                device_scratch_padA, device_velocity, device_position,
+                array_size, &device_ray_direction[i*3]);
+
+    }
+
+    CHECK_CUDA_ERROR(cudaFree(dev_dot_product_result));
+    CHECK_CUDA_ERROR(cudaFree(dev_direction_vec));
+
+}
+
+
+void universe::process_collisions_with_other_particles(void)
+{
+    
+}
+
+void universe::process_particles(void)
+{
+    const float delta = DELTA_TIME;
+    
+    // // Update the positions of the particles:
+    // // self.positions[counti] = self.velocities[counti] * DELTA_TIME + self.positions[counti]
+    CHECK_CUDA_ERROR(cublasSaxpy(cublas_handle, array_size*3, &delta,
+                device_velocity, 1, device_position, 1));
+
+    // Update the velocities of the particles:
+    // self.velocities[counti] = self.accelerations[counti] * DELTA_TIME + self.velocities[counti]
+    CHECK_CUDA_ERROR(cublasSaxpy(cublas_handle, array_size*3, &delta,
+                device_acceleration, 1, device_velocity, 1));
+
+    process_collisions_with_rays();
+
+    process_collisions_with_other_particles();
+}
+
+void universe::get_positions(float ** in_out_buffer, int * size, int * stride)
+{
+    CHECK_CUDA_ERROR(cudaMemcpy(host_position, device_position, 
+                long(array_size * 3 * sizeof(float)), cudaMemcpyDeviceToHost));
+
+    *in_out_buffer = host_position;
+    *size = num_of_particles;
+    *stride = array_size;
+}
+
+universe::universe(int size)
+{
+curandCreateGenerator(&gen,
+            CURAND_RNG_PSEUDO_DEFAULT);
+
+    CHECK_CUDA_ERROR(cublasCreate(&cublas_handle));
+
+    num_of_particles = size;
+    array_size = size; 
+    
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(device_position), 
+                long(array_size*3*sizeof(float))));
+
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(device_velocity), 
+                long(array_size*3*sizeof(float))));
+
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(device_acceleration),
+                long(array_size*3*sizeof(float))));
+
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(device_scratch_padA),
+                long(array_size*3*sizeof(float))));
+
+    host_position = (float *)malloc(long(array_size*3*sizeof(float)));
+    host_velocity = (float *)malloc(long(array_size*3*sizeof(float)));
+    host_acceleration = (float *)malloc(long(array_size*3*sizeof(float)));
+
+    CHECK_CUDA_ERROR(curandGenerateUniform(gen, device_position, array_size*2));
+    CHECK_CUDA_ERROR(nppsAddC_32f_I(20.0, device_position, array_size*3));
+    CHECK_CUDA_ERROR(curandGenerateUniform(gen, device_velocity, array_size*2));
+    // CHECK_CUDA_ERROR(curandGenerateUniform(gen, device_acceleration, num_of_particles*3));
+
+    //Add rays that can collied with particles to the universe
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(device_ray_position),
+                long(4*3*sizeof(float))));
+
+    CHECK_CUDA_ERROR(cudaMallocManaged(&(device_ray_direction),
+                long(4*3*sizeof(float))));
+    
+    float temp_positions[] = {0.0,0.0,0.0,0.0,0.0,0.0,800.0,800.0,800.0,800.0,800.0,800.0};
+    float temp_directions[] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0};
+
+    CHECK_CUDA_ERROR(cudaMemcpy(device_ray_position, temp_positions,
+                long(sizeof(temp_positions)), cudaMemcpyHostToDevice));
+
+    CHECK_CUDA_ERROR(cudaMemcpy(device_ray_direction, temp_directions,
+                long(sizeof(temp_directions)), cudaMemcpyHostToDevice));
+
+    
+}
+
+universe::~universe()
+{
+
+
+    CHECK_CUDA_ERROR(cudaFree(device_position));
+    CHECK_CUDA_ERROR(cudaFree(device_velocity));
+    CHECK_CUDA_ERROR(cudaFree(device_acceleration));
+    CHECK_CUDA_ERROR(cudaFree(device_scratch_padA));
+
+    free(host_position);
+    free(host_velocity);
+    free(host_acceleration);
+    CHECK_CUDA_ERROR(cublasDestroy(cublas_handle));
+
+
+
+}
+
+void universe::process(void)
+{
+
+    //First we need to process gravity
+
+
+
+    //Then we want to process collisions
+    process_particles();
+
+}

Cpp_Stuff/src/utility/custom_kernels.cu

@@ -0,0 +1,39 @@
+#include "custom_kernels.cuh"
+#include "stdio.h"
+#define BLOCK_SIZE 32
+
+__global__ static void process_ray_collision(int size, float * norms, float * direction, float * vel, float * pos, int array_stride, float * line_vec)
+{
+    int index = blockIdx.x * blockDim.x + threadIdx.x;
+    int stride = blockDim.x * gridDim.x;
+    for (int i = index; i < size; i+= stride)
+    {
+        if(norms[i] < 5.0)
+        {
+            float speed_norm = 0.0;
+            float dot_product = 0.0;
+            for(int k = 0; k<3; k++)
+            {
+                dot_product += vel[k*array_stride+i]*line_vec[k];
+            }
+            for(int k = 0; k<3; k++)
+            {
+                float temp = (dot_product*(line_vec[k]));
+                speed_norm += pow(temp, 2); 
+            } 
+            speed_norm = sqrt(speed_norm)*2.0;
+            for(int k = 0; k<3; k++)
+            {
+                vel[k*array_stride+i] = vel[k*array_stride+i] - speed_norm*direction[k*size+i];
+                pos[k*array_stride+i] = pos[k*array_stride+i] - direction[k*size+i]*(5.0-norms[i]);
+            }
+        }
+    }
+}
+
+
+void custom_kernel_process_ray_collision(int size, float * norms, float * direction, float * vel, float * pos, int array_stride, float * line_vec)
+{
+    int numBlocks = (size + BLOCK_SIZE - 1) / BLOCK_SIZE;
+    process_ray_collision<<<numBlocks,BLOCK_SIZE>>>(size, norms, direction, vel, pos, array_stride, line_vec);
+}

Cpp_Stuff/src/utility/hal.cpp

@@ -0,0 +1,25 @@
+#include "hal.hpp"
+#include "npp.h"
+#include <cuda.h>
+#include <cuda_runtime.h>
+
+std::vector<float> hal_vector_f32_cross_multiply(std::vector<float>& a, std::vector<float>& b)
+{
+    Npp32f * pSrca; 
+    Npp32f * pSrcb; 
+    cudaMallocManaged((void **)(&pSrca), long(sizeof(Npp32f) * a.size()));
+    cudaMallocManaged((void **)(&pSrcb), long(sizeof(Npp32f) * b.size()));
+    cudaMemcpy(pSrca, &a[0], sizeof(Npp32f) * a.size(), cudaMemcpyHostToDevice);
+    cudaMemcpy(pSrcb, &b[0], sizeof(Npp32f) * b.size(), cudaMemcpyHostToDevice);
+    std::vector<float> ret_vec;
+    ret_vec = b;
+    Npp32f * pSrcc; 
+    cudaMallocManaged((void **)(&pSrcc), long(sizeof(Npp32f) * ret_vec.size()));
+    NppStatus status = nppsMul_32f(pSrca, pSrcb, pSrcc, a.size());
+    cudaDeviceSynchronize();
+    cudaMemcpy(&ret_vec[0], pSrcc, sizeof(Npp32f) * ret_vec.size(), cudaMemcpyDeviceToHost);
+    cudaFree(pSrca);
+    cudaFree(pSrcb);
+    cudaFree(pSrcc);
+    return ret_vec;
+}

Python_Stuff/python_particle_simulator.py

@@ -0,0 +1,62 @@
+import numpy as np
+import time
+import sys
+import sdl2
+import sdl2.ext
+
+import world_objects
+import universe
+import renderer
+
+sprites = ["./sdl_sprites/green.bmp",
+"./sdl_sprites/red.bmp",
+"./sdl_sprites/yellow.bmp",
+"./sdl_sprites/blue.bmp"]
+        
+def generate_particle(my_renderer):
+    location = np.random.randint(10,750,2)
+    accel = np.random.randint(0,5,2)
+    vel = np.random.randint(-10,10,2)
+    print(len(sprites))
+    sprite = sprites[np.random.randint(0,(len(sprites)))]
+    return world_objects.particle(location, vel, accel, 10, (255,0,0), 10.0, "blah", my_renderer, sprite)
+    
+        
+def run():
+    my_renderer = renderer.SdlRender()
+    running = True
+    my_universe = universe.universe()
+    black_hole = world_objects.particle([400,400], [10,0], [0,0], 20, (0,0,255), 2000.0, "red", my_renderer, sprites[1])
+    my_universe.add_particle(black_hole)
+
+    for i in range(100):
+        my_universe.add_particle(generate_particle(my_renderer))
+    
+    left = world_objects.line([0,0], [1,0])
+    top = world_objects.line([0,0], [0,1])
+    right = world_objects.line([800, 800], [1,0])
+    bottom = world_objects.line([800, 800], [0,1])
+    my_universe.add_line(left)
+    my_universe.add_line(top)
+    my_universe.add_line(right)
+    my_universe.add_line(bottom)
+    
+    print(my_universe)
+    start_time=0
+    while running:
+        start_time = time.time()
+        events = sdl2.ext.get_events()
+        for event in events:
+            if event.type == sdl2.SDL_QUIT:
+                running = False
+                break
+        my_universe.process()
+        stop_time = time.time()
+        
+        my_renderer.sdl_world.process()
+        temp = stop_time - start_time
+        start_time = time.time()
+        print(f"Process physics: {(temp)*1000} ms Draw world: {(start_time - stop_time)*1000} ms")
+        
+if __name__ == '__main__':
+    sys.exit(run())

Python_Stuff/renderer.py

@@ -0,0 +1,37 @@
+import sdl2
+import sdl2.ext
+from sdl2.ext.image import load_image
+import numpy as np
+
+
+WHITE = sdl2.ext.Color(255, 255, 255)
+
+class Player(sdl2.ext.Entity):
+    def __init__(self, world, renderer, radius, path_to_image, posx=0, posy=0):
+        image = load_image(str(path_to_image))
+        new_surf = sdl2.surface.SDL_CreateRGBSurface(0,radius*2,radius*2,32,0x000000FF, 0x0000FF00, 0x00FF0000, 0xFF000000)
+        sdl2.surface.SDL_BlitScaled(image, sdl2.rect.SDL_Rect(0,0,image.h,image.w), new_surf, sdl2.rect.SDL_Rect(0,0,radius*2,radius*2))
+        self.sprite = renderer.factory.from_surface(new_surf, False)
+        self.sprite.position = int(posx), int(posy)
+
+    def update_pos(self, pos):
+        self.sprite.position = int(pos[0]), int(pos[1])
+        
+class SoftwareRenderer(sdl2.ext.SoftwareSpriteRenderSystem):
+    def __init__(self, window):
+        super(SoftwareRenderer, self).__init__(window)
+
+    def render(self, components):
+        sdl2.ext.fill(self.surface, sdl2.ext.Color(0, 0, 0))
+        super(SoftwareRenderer, self).render(components)
+        
+class SdlRender():
+    def __init__(self):
+        sdl2.ext.init()
+        self.window = sdl2.ext.Window("Gravity Sim", size=(800, 800))
+        self.window.show()
+        self.sdl_world = sdl2.ext.World()
+
+        self.spriterenderer = SoftwareRenderer(self.window)
+        self.sdl_world.add_system(self.spriterenderer)
+        self.factory = sdl2.ext.SpriteFactory(sdl2.ext.SOFTWARE)

Python_Stuff/requirements.txt

@@ -0,0 +1,2 @@
+numpy==1.24.1
+PySDL2==0.9.14

Python_Stuff/sdl.py

@@ -0,0 +1,50 @@
+import sys
+import sdl2
+import sdl2.ext
+
+WHITE = sdl2.ext.Color(255, 255, 255)
+
+class SoftwareRenderer(sdl2.ext.SoftwareSpriteRenderSystem):
+    def __init__(self, window):
+        super(SoftwareRenderer, self).__init__(window)
+
+    def render(self, components):
+        sdl2.ext.fill(self.surface, sdl2.ext.Color(0, 0, 0))
+        super(SoftwareRenderer, self).render(components)
+
+
+class Player(sdl2.ext.Entity):
+    def __init__(self, world, sprite, posx=0, posy=0):
+        self.sprite = sprite
+        self.sprite.position = posx, posy
+        
+def run():
+    sdl2.ext.init()
+    window = sdl2.ext.Window("The Pong Game", size=(800, 600))
+    window.show()
+    world = sdl2.ext.World()
+
+    spriterenderer = SoftwareRenderer(window)
+    world.add_system(spriterenderer)
+    factory = sdl2.ext.SpriteFactory(sdl2.ext.SOFTWARE)
+
+    sp_paddle1 = factory.from_color(WHITE, size=(20, 100))
+    sp_paddle2 = factory.from_color(WHITE, size=(20, 100))
+
+    player1 = Player(world, sp_paddle1, 0, 250)
+    player2 = Player(world, sp_paddle2, 780, 250)
+
+    running = True
+    while running:
+        events = sdl2.ext.get_events()
+        for event in events:
+            if event.type == sdl2.SDL_QUIT:
+                running = False
+                break
+        
+        world.process()
+        #window.refresh()
+    return 0
+
+if __name__ == "__main__":
+    sys.exit(run())

Python_Stuff/universe.py

@@ -0,0 +1,153 @@
+
+import numpy as np
+
+DELTA_TIME = 0.05
+GRAV_CONST = 6.674  # m^3*s^2/kg
+
+
+class universe():
+    def __init__(self):
+        self.matA = None
+        self.matB = None
+        self.matC = None
+        self.line_array = []
+        self.mass_array = []
+        self.radius_array = []
+        self.force_direction = []
+        self.positions = []
+        self.velocities = []
+        self.accelerations = []
+        self.particle_list = []
+        self.sprites = []
+
+        print("New universe has been created")
+
+    def __repr__(self):
+        return "Your universe contains the following particles: " +str(self.particle_list)
+
+    def _recalculate_matA(self):
+        #
+        #      | 0   m2   m3   m4  |
+        #      |                               |
+        #      | m1  0    m3   m4  |
+        #  A = |                               | * G
+        #      | m1  m2   0    m4  |
+        #      |                               |
+        #      | m1  m2   m3   0   |
+        #      |                               |
+        #
+        temp_array = np.matrix(self.mass_array)
+        self.matA = np.matmul(np.ones(temp_array.T.shape), temp_array)
+        self.matA = GRAV_CONST * self.matA
+        self.matA = self.matA - np.multiply(self.matA, np.eye(*self.matA.shape))
+        print(self.matA)
+
+    def _recalculate_matB(self):
+        #      | 1          [x1-x2]^2   [x1-x3]^2   [x1-x4]^2  |
+        #      |                                               |
+        #      | [x2-x1]^2  1           [x2-x3]^2   [x2-x4]^2  |
+        #  B = |                                               |
+        #      | [x3-x1]^2  [x3-x2]^2   1           [x3-x4]^2  |
+        #      |                                               |
+        #      | [x4-x1]^2  [x4-x2]^2   [x4-x3]^2   1          |
+        #      |                                               |
+
+        temp_array = np.tile(np.array(self.positions), (len(self.positions), 1, 1))
+        self.force_direction = temp_array - np.transpose(temp_array, axes=(1, 0, 2))
+        self.matB = np.empty(self.force_direction.shape[0:2])
+        for counti, i in enumerate(self.force_direction):
+            for countj, j in enumerate(i):
+                self.matB[counti, countj] = np.dot(j, j)
+                if self.matB[counti, countj] == 0.0:
+                    self.matB[counti, countj] = 1.0
+                self.force_direction[counti, countj] = j / np.sqrt(self.matB[counti, countj])
+
+    def _recalculate_matC(self):
+        #      | 0          [r1+r2]^2   [r1+r3]^2   [r1+r4]^2  |
+        #      |                                               |
+        #      | [r2+r1]^2  0           [r2+r3]^2   [r2+r4]^2  |
+        #  C = |                                               |
+        #      | [r3+r1]^2  [r3+r2]^2   0           [r3+r4]^2  |
+        #      |                                               |
+        #      | [r4+r1]^2  [r4+r2]^2   [r4+r3]^2   0          |
+        #      |                                               |
+
+        temp_array = np.matrix(self.radius_array)
+        self.matC = np.matmul(np.ones(temp_array.T.shape), temp_array) + temp_array.T
+        self.matC = np.multiply(self.matC, self.matC)
+        # print("Matrix C:", self.matC)
+
+    def add_particle(self, particle):
+        self.mass_array.append(particle.mass)
+        self.positions.append(particle.position)
+        self.velocities.append(particle.velocity)
+        self.accelerations.append(particle.acceleration)
+        self.particle_list.append(particle)
+        self.radius_array.append(particle.radius)
+        self.sprites.append(particle.sprite)
+        # print("Added new particle")
+        self._recalculate_matA()
+        self._recalculate_matB()
+        self._recalculate_matC()
+
+    def add_line(self, line):
+        self.line_array.append(line)
+
+    # def process_collisions(self):
+
+    # def update_gravity(self):
+
+    def process(self):
+
+        # Apply gravity
+        temp = self.matA/self.matB
+        for counti in range(temp.shape[0]):
+            accel = np.array([0.0, 10.0])
+            for countj in range(temp.shape[1]):
+                accel += temp[counti, countj] * self.force_direction[counti, countj]
+
+            self.accelerations[counti] = accel
+            self.velocities[counti] = self.accelerations[counti] * DELTA_TIME + self.velocities[counti]
+            self.positions[counti] = self.velocities[counti] * DELTA_TIME + self.positions[counti]
+            self.sprites[counti].update_pos(self.positions[counti]-self.radius_array[counti])
+        self._recalculate_matB()
+
+        # Check for collisions...
+        for counti in range(temp.shape[0]):
+            # ...with other particles
+            for countj in range(counti+1, temp.shape[0]):
+                if (self.matB[counti, countj] <= self.matC[counti, countj]):
+                    direction_vector = self.positions[countj] - self.positions[counti]
+                    dist = np.linalg.norm(direction_vector)
+                    direction_vector = direction_vector / dist
+                    temp_vel_1 = np.dot(self.velocities[counti], direction_vector)*(direction_vector)
+                    leftover_vel_1 = self.velocities[counti] - temp_vel_1
+                    temp_vel_2 = np.dot(self.velocities[countj], direction_vector)*(direction_vector)
+                    leftover_vel_2 = self.velocities[countj] - temp_vel_2
+                    m1_2 = self.mass_array[counti]*2  
+                    m2_2 = self.mass_array[countj]*2  
+                    m1_plus_m2 = self.mass_array[counti]+self.mass_array[countj]
+                    m1_minus_m2 = self.mass_array[counti] - self.mass_array[countj]
+                    self.velocities[countj] = (0.9999*m1_2*temp_vel_1-m1_minus_m2*temp_vel_2)/m1_plus_m2 + leftover_vel_2
+                    self.velocities[counti] = (0.9999*m1_minus_m2*temp_vel_1+m2_2*temp_vel_2)/m1_plus_m2 + leftover_vel_1
+
+                    # self.positions[countj] = self.velocities[countj] * DELTA_TIME + self.positions[countj]
+                    # self.positions[counti] = self.velocities[counti] * DELTA_TIME + self.positions[counti]
+
+                    self.positions[countj] = self.positions[countj] + direction_vector*((np.sqrt(self.matC[counti,countj])-np.sqrt(self.matB[counti,countj]))*self.mass_array[counti]/m1_plus_m2)
+                    self.positions[counti] = self.positions[counti] - direction_vector*((np.sqrt(self.matC[counti,countj])-np.sqrt(self.matB[counti,countj]))*self.mass_array[countj]/m1_plus_m2)
+
+            # with lines
+            for line in self.line_array:
+                vec1 = line.position - np.array(self.positions[counti])
+                dist = np.dot(vec1, line.vector)*line.vector
+                norm = np.linalg.norm(dist)
+                direction = dist/norm
+                if(norm < self.radius_array[counti]):
+                    # print("HIT!")
+                    speed = np.dot(self.velocities[counti], line.vector)*line.vector
+                    self.velocities[counti] = self.velocities[counti] - 2*np.linalg.norm(speed)*direction
+
+                    # self.positions[counti] = self.velocities[counti] * DELTA_TIME + self.positions[counti]
+                    self.positions[counti] = self.positions[counti]-direction*(self.radius_array[counti]-norm)
+

Python_Stuff/world_objects.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+from renderer import Player
+
+
+class particle():
+    def __init__(self, pos, vel, accel, radius, color, mass, name, renderer, path_to_image):
+        self.position = np.array(pos, dtype=np.float32)
+        self.acceleration = np.array(accel, dtype=np.float32)
+        self.velocity = np.array(vel, dtype=np.float32)
+        self.radius = radius
+        self.color = color
+        self.mass = mass
+        self.name = name
+        self.sprite = Player(renderer.sdl_world, renderer, radius, path_to_image, pos[0], pos[1])
+        
+class line():
+    def __init__(self, pos, vector):
+        self.position = np.array(pos, dtype=np.float32)
+        self.vector = np.array(vector, dtype=np.float32)