Raytracer slice in cuda, end and start angles as output

Co-Authored-By: Emil Wallberg <49481622+EmilWallberg@users.noreply.github.com>

modules/blackhole/cuda/blackhole_cuda.cu

@@ -2,76 +2,89 @@
 #include <vector>
 #include "device_launch_parameters.h"
 
-__constant__ double rs = 1;
+__constant__ double PI = 3.1415926535897932384626433832795;
 
-__device__ void geodesic(double u, double du, double& out_u, double& out_du) {
-    double temp_u = u;
-    out_u = du;
-    out_du = -temp_u * (1 - (3 / 2) * rs * temp_u);
+__device__ void geodesic(double u, double dudphi, double& out_du_dphi, double& out_d2u_dphi2, double rs) {
+    out_du_dphi = dudphi;
+    out_d2u_dphi2 = -u * (1 - 3. / 2. * rs * u);
 }
 
-__device__ void rk4_step(double& u, double& du, double& phi, double h) {
-    double k1_u, k1_du, k2_u, k2_du, k3_u, k3_du, k4_u, k4_du;
+__device__ void rk4_step(double& u, double& dudphi, double& phi, double h, double rs) {
+    double k1_u, k1_dudphi, k2_u, k2_dudphi, k3_u, k3_dudphi, k4_u, k4_dudphi;
 
-    geodesic(u, du, k1_u, k1_du);
-    geodesic(u + 0.5 * k1_u * h, du + 0.5 * k1_du * h, k2_u, k2_du);
-    geodesic(u + 0.5 * k2_u * h, du + 0.5 * k2_du * h, k3_u, k3_du);
-    geodesic(u + k3_u * h, du + k3_du * h, k4_u, k4_du);
+    geodesic(u, dudphi, k1_u, k1_dudphi, rs);
+    geodesic(u + 0.5 * k1_u * h, dudphi + 0.5 * k1_dudphi * h, k2_u, k2_dudphi, rs);
+    geodesic(u + 0.5 * k2_u * h, dudphi + 0.5 * k2_dudphi * h, k3_u, k3_dudphi, rs);
+    geodesic(u + k3_u * h, dudphi + k3_dudphi * h, k4_u, k4_dudphi, rs);
 
     phi += h;
     u = u + (k1_u + 2 * k2_u + 2 * k3_u + k4_u) * h / 6;
-    du = du + (k1_du + 2 * k2_du + 2 * k3_du + k4_du) * h / 6;
+    dudphi = dudphi + (k1_dudphi + 2 * k2_dudphi + 2 * k3_dudphi + k4_dudphi) * h / 6;
 }
 
-__global__ void solve_geodesic_kernel(double u_0, double* du_0_values, double h, int num_paths, int num_steps, double* u_values, double* du_values, double* phi_values) {
+__global__ void solveGeodesicKernel(double rs, double envmap_r, double u_0, double* dudphi_0_values, double h, int num_paths, int num_steps, double* angles_out) {
     int idx = blockIdx.x * blockDim.x + threadIdx.x;
     if (idx >= num_paths) return;
 
     double u = u_0;
-    double du = du_0_values[idx];
+    double dudphi = dudphi_0_values[idx];
     double phi = 0.0;
 
-    for (int step = 0; step < num_steps; step++) {
-        u_values[idx * num_steps + step] = u;
-        du_values[idx * num_steps + step] = du;
-        phi_values[idx * num_steps + step] = phi;
+    // Perform the first RK4 step
+    rk4_step(u, dudphi, phi, h, rs);
 
-        rk4_step(u, du, phi, h);
+    double r = 1.0 / u;
+    double r_0 = 1.0 / u_0;
+    double a = r * sin(phi);
+    double b = r_0 - r * cos(phi);
+
+    // Store starting angle (angle is local to camera)
+    angles_out[idx * 2] = atan2(a, b);
+
+    auto out_of_bounds = [&u, &envmap_r]() -> bool {
+        return (1.0 / u > envmap_r);
+        };
+
+    auto inside_singularity = [&u, &rs]() -> bool {
+        return (1.0 / u <= rs);
+        };
+
+    for (int step = 1; step < num_steps && !out_of_bounds() && !inside_singularity(); step++) {
+        rk4_step(u, dudphi, phi, h, rs);
+    }
+
+    angles_out[idx * 2 + 1] = !inside_singularity() ? phi : nan("");
+}
+
+void generate_du(double* d_du_0_values, double min, double max, int count) {
+    double step = (max - min) / (count - 1);
+    for (int i = 0; i < count; ++i) {
+        d_du_0_values[i] = min + step * i;
     }
 }
 
-extern "C" {
-    __declspec(dllexport) void cuda_test(
-        int num_paths, int num_steps, double u_0,
-        double* du_0_values, double h, double* u_out, double* phi_out) {
+void schwarzchild(
+    double rs, double envmap_r, int num_paths, int num_steps, double u_0, double h, double* angle_out) {
 
-        // Allocate device memory
-        double* d_du_0_values;
-        double* d_u_values;
-        double* d_du_values;
-        double* d_phi_values;
+    double* d_dudphi_0_values;
+    double* d_angle_values;
 
-        cudaMalloc(&d_du_0_values, num_paths * sizeof(double));
-        cudaMalloc(&d_u_values, num_paths * num_steps * sizeof(double));
-        cudaMalloc(&d_du_values, num_paths * num_steps * sizeof(double));
-        cudaMalloc(&d_phi_values, num_paths * num_steps * sizeof(double));
+    // Allocate device memory
+    cudaMalloc(&d_dudphi_0_values, num_paths * sizeof(double));
+    cudaMalloc(&d_angle_values, num_paths * 2 * sizeof(double));
 
+    // Copy initial velocity values to device
+    std::vector<double> dudphi_0_values(num_paths, 0);
+    generate_du(dudphi_0_values.data(), 1, -1, num_paths);
+    cudaMemcpy(d_dudphi_0_values, dudphi_0_values.data(), num_paths * sizeof(double), cudaMemcpyHostToDevice);
 
-        // Copy initial velocity values to device
-        cudaMemcpy(d_du_0_values, du_0_values, num_paths * sizeof(double), cudaMemcpyHostToDevice);
+    // Launch kernel
+    int threadsPerBlock = 256;
+    int numBlocks = (num_paths + threadsPerBlock - 1) / threadsPerBlock;
+    solveGeodesicKernel<<<numBlocks, threadsPerBlock>>>(rs, envmap_r, u_0, d_dudphi_0_values, h, num_paths, num_steps, d_angle_values);
 
-        // Launch kernel
-        int threadsPerBlock = 256;
-        int numBlocks = (num_paths + threadsPerBlock - 1) / threadsPerBlock;
-        solve_geodesic_kernel << <numBlocks, threadsPerBlock >> > (u_0, d_du_0_values, h, num_paths, num_steps, d_u_values, d_du_values, d_phi_values);
+    cudaMemcpy(angle_out, d_angle_values, num_paths * 2 * sizeof(double), cudaMemcpyDeviceToHost);
 
-        // Copy results back to host
-        cudaMemcpy(u_out, d_u_values, num_paths * num_steps * sizeof(double), cudaMemcpyDeviceToHost);
-        cudaMemcpy(phi_out, d_phi_values, num_paths * num_steps * sizeof(double), cudaMemcpyDeviceToHost);
-
-        cudaFree(d_du_0_values);
-        cudaFree(d_u_values);
-        cudaFree(d_du_values);
-        cudaFree(d_phi_values);
-    }
+    cudaFree(d_dudphi_0_values);
+    cudaFree(d_angle_values);
 }

modules/blackhole/cuda/blackhole_cuda.h

@@ -2,7 +2,6 @@
 #define __OPENSPACE_MODULE_BLACKHOLE___CUDA___H__
 #include <vector>
 
-void cuda_test(int num_paths, int num_steps, double u_0, std::vector<double>& du_0_values, std::vector<double>& u_out, std::vector<double>& phi_out);
-
+void schwarzchild(double rs, double envmap_r, int num_paths, int num_steps, double u_0, double h, double* angle_out);
 
 #endif

modules/blackhole/rendering/renderableblackhole.cpp

@@ -12,6 +12,7 @@
 #include <ghoul/logging/logmanager.h>
 #include <ghoul/filesystem/filesystem.h>
 #include <filesystem>
+#include <vector>
 
 #include <modules/blackhole/cuda/blackhole_cuda.h>
 
@@ -37,7 +38,11 @@ namespace openspace {
     void RenderableBlackHole::initialize() {
         _viewport = ViewPort(global::navigationHandler->camera());
         global::navigationHandler->camera()->setRotation(glm::dquat(0,0,0,0));
-        //cuda_test();
+        _schwarzschildWarpTable = std::vector<double>(_rayCount * 2, 0);
+        schwarzchild(1, 30, _rayCount, 5000, 1.0 / 20.0, 0.01, _schwarzschildWarpTable.data());
+        for (int i = 0; i < _rayCount; ++i) {
+            LINFO(std::format("phi = [s: {:.10f}, e: {:.10f}]", _schwarzschildWarpTable[i * 2], _schwarzschildWarpTable[i * 2 + 1]));
+        }
     }
 
     void RenderableBlackHole::initializeGL() {
@@ -122,8 +127,10 @@ namespace openspace {
     }
 
     void RenderableBlackHole::loadEnvironmentTexture() {
-        const std::string texturePath = "${MODULE_BLACKHOLE}/rendering/uv.png";
+        //const std::string texturePath = "${MODULE_BLACKHOLE}/rendering/uv.png";
         //const std::string texturePath = "C:/Users/wilbj602/Documents/GitHub/OpenSpace/sync/http/milkyway_textures/2/DarkUniverse_mellinger_8k.jpg";
+        const std::string texturePath = "C:/Users/wilbj602/Downloads/img.jpg";
+
         _enviromentTexture = ghoul::io::TextureReader::ref().loadTexture(absPath(texturePath), 2);
 
         if (_enviromentTexture) {

modules/blackhole/rendering/renderableblackhole.h

@@ -34,7 +34,8 @@ private:
     void loadEnvironmentTexture();
   ghoul::opengl::ProgramObject* _program = nullptr;
   ViewPort _viewport;
-
+  size_t _rayCount = 100;
+  std::vector<double> _schwarzschildWarpTable;
 
   GLuint _framebuffer = 0;
   GLuint _quadVao = 0;