From 6eaa87c433bbbab3cfb7b644935071b614e28975 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Wilhelm=20Bj=C3=B6rkstr=C3=B6m?= <143391787+Grantallkotten@users.noreply.github.com> Date: Wed, 5 Mar 2025 11:22:42 +0100 Subject: [PATCH] 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 | 107 ++++++++++-------- modules/blackhole/cuda/blackhole_cuda.h | 3 +- .../rendering/renderableblackhole.cpp | 11 +- .../blackhole/rendering/renderableblackhole.h | 3 +- 4 files changed, 72 insertions(+), 52 deletions(-) diff --git a/modules/blackhole/cuda/blackhole_cuda.cu b/modules/blackhole/cuda/blackhole_cuda.cu index ba994718c0..68f82c88ab 100644 --- a/modules/blackhole/cuda/blackhole_cuda.cu +++ b/modules/blackhole/cuda/blackhole_cuda.cu @@ -2,76 +2,89 @@ #include #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 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<<>>(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 << > > (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); } diff --git a/modules/blackhole/cuda/blackhole_cuda.h b/modules/blackhole/cuda/blackhole_cuda.h index 2ab200bdda..f2314e93c1 100644 --- a/modules/blackhole/cuda/blackhole_cuda.h +++ b/modules/blackhole/cuda/blackhole_cuda.h @@ -2,7 +2,6 @@ #define __OPENSPACE_MODULE_BLACKHOLE___CUDA___H__ #include -void cuda_test(int num_paths, int num_steps, double u_0, std::vector& du_0_values, std::vector& u_out, std::vector& phi_out); - +void schwarzchild(double rs, double envmap_r, int num_paths, int num_steps, double u_0, double h, double* angle_out); #endif diff --git a/modules/blackhole/rendering/renderableblackhole.cpp b/modules/blackhole/rendering/renderableblackhole.cpp index 2b2ac7e206..2602eeb48f 100644 --- a/modules/blackhole/rendering/renderableblackhole.cpp +++ b/modules/blackhole/rendering/renderableblackhole.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include @@ -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(_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) { diff --git a/modules/blackhole/rendering/renderableblackhole.h b/modules/blackhole/rendering/renderableblackhole.h index 81a3be9b11..aa367d26e7 100644 --- a/modules/blackhole/rendering/renderableblackhole.h +++ b/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 _schwarzschildWarpTable; GLuint _framebuffer = 0; GLuint _quadVao = 0;