diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..9b6e850 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +build +.vscode +data \ No newline at end of file diff --git a/makefile b/Makefile similarity index 74% rename from makefile rename to Makefile index 18993ff..2123394 100644 --- a/makefile +++ b/Makefile @@ -1,6 +1,6 @@ # Compiler and flags NVCC = nvcc -CXXFLAGS = -I./src -I./linalg -I./img -I./objs -std=c++17 +CXXFLAGS = -I./src -I./hurricanedata -std=c++17 $(shell nc-config --cxx4flags) $(shell nc-config --cxx4libs) # Directories SRC_DIR = src @@ -8,7 +8,7 @@ BUILD_DIR = build # Files TARGET = $(BUILD_DIR)/main -SRC_FILES = $(wildcard $(SRC_DIR)/*.cu) +SRC_FILES := $(shell find $(SRC_DIR) -type f \( -name '*.cu' -o -name '*.cpp' \)) OBJ_FILES = $(patsubst $(SRC_DIR)/%.cu,$(BUILD_DIR)/%.o,$(SRC_FILES)) # Default target @@ -20,6 +20,7 @@ $(TARGET): $(OBJ_FILES) | $(BUILD_DIR) # Compile object files $(BUILD_DIR)/%.o: $(SRC_DIR)/%.cu + @mkdir -p $(dir $@) $(NVCC) $(CXXFLAGS) -c $< -o $@ # Debug build diff --git a/build/main b/build/main index ffdc0f6..4093b44 100755 Binary files a/build/main and b/build/main differ diff --git a/build/main.o b/build/main.o index c0cf63c..d5059a5 100644 Binary files a/build/main.o and b/build/main.o differ diff --git a/output.ppm b/output.ppm deleted file mode 100644 index 1d7b780..0000000 Binary files a/output.ppm and /dev/null differ diff --git a/src/hurricanedata/datareader.cu b/src/hurricanedata/datareader.cu new file mode 100644 index 0000000..6052937 --- /dev/null +++ b/src/hurricanedata/datareader.cu @@ -0,0 +1,25 @@ +#include "datareader.h" + +#include + +using namespace std; +using namespace netCDF; + +std::vector readData(std::string path, std::string variableName) { + netCDF::NcFile data(path, netCDF::NcFile::read); + + multimap vars = data.getVars(); + + NcVar var = vars.find(variableName)->second; + + int length = 1; + for (NcDim dim: var.getDims()) { + length *= dim.getSize(); + } + + vector vec(length); + + var.getVar(vec.data()); + + return vec; +} \ No newline at end of file diff --git a/src/hurricanedata/datareader.h b/src/hurricanedata/datareader.h new file mode 100644 index 0000000..7e628cc --- /dev/null +++ b/src/hurricanedata/datareader.h @@ -0,0 +1,9 @@ +#ifndef DATAREADER_H +#define DATAREADER_H + +#include +#include + +std::vector readData(std::string path, std::string variableName); + +#endif //DATAREADER_H diff --git a/src/img/handler.h b/src/img/handler.h deleted file mode 100644 index 9fe80f5..0000000 --- a/src/img/handler.h +++ /dev/null @@ -1,13 +0,0 @@ -#pragma once - -#include - - -void saveImage(const char* filename, unsigned char* framebuffer, int width, int height) { // TODO: Figure out a better way to do this - std::ofstream imageFile(filename, std::ios::out | std::ios::binary); - imageFile << "P6\n" << width << " " << height << "\n255\n"; - for (int i = 0; i < width * height * 3; i++) { - imageFile << framebuffer[i]; - } - imageFile.close(); -} \ No newline at end of file diff --git a/src/linalg/linalg.h b/src/linalg/linalg.h deleted file mode 100644 index 885273e..0000000 --- a/src/linalg/linalg.h +++ /dev/null @@ -1,9 +0,0 @@ -#pragma once - -#ifndef LINALG_H -#define LINALG_H - -#include "vec.h" -#include "mat.h" - -#endif // LINALG_H diff --git a/src/linalg/mat.h b/src/linalg/mat.h deleted file mode 100644 index 6f70f09..0000000 --- a/src/linalg/mat.h +++ /dev/null @@ -1 +0,0 @@ -#pragma once diff --git a/src/linalg/vec.h b/src/linalg/vec.h deleted file mode 100644 index ef505bd..0000000 --- a/src/linalg/vec.h +++ /dev/null @@ -1,18 +0,0 @@ -#pragma once - -#include -#include - -struct Vec3 { // TODO: Maybe make this into a class - double x, y, z; - - __host__ __device__ Vec3() : x(0), y(0), z(0) {} - __host__ __device__ Vec3(double x, double y, double z) : x(x), y(y), z(z) {} - - __host__ __device__ Vec3 operator+(const Vec3& b) const { return Vec3(x + b.x, y + b.y, z + b.z); } - __host__ __device__ Vec3 operator-(const Vec3& b) const { return Vec3(x - b.x, y - b.y, z - b.z); } - __host__ __device__ Vec3 operator*(double b) const { return Vec3(x * b, y * b, z * b); } - __host__ __device__ Vec3 operator-() const { return Vec3(-x, -y, -z); } - __host__ __device__ double dot(const Vec3& b) const { return x * b.x + y * b.y + z * b.z; } - __host__ __device__ Vec3 normalize() const { double len = sqrt(x * x + y * y + z * z); return Vec3(x / len, y / len, z / len); } -}; \ No newline at end of file diff --git a/src/main.cu b/src/main.cu index fafaec8..ccd7aef 100644 --- a/src/main.cu +++ b/src/main.cu @@ -1,102 +1,21 @@ +#include "hurricanedata/datareader.h" + #include #include #include #include -#include "linalg/linalg.h" -#include "objs/sphere.h" -#include "img/handler.h" +int main() { + std::string path = "data/MERRA2_400.inst6_3d_ana_Np.20120101.nc4"; + std::string variable = "U"; + auto x = readData(path, variable); -#define WIDTH 3840 -#define HEIGHT 2160 -#define SAMPLES_PER_PIXEL 8 - - -__device__ Vec3 phongShading(const Vec3& point, const Vec3& normal, const Vec3& lightDir, const Vec3& viewDir, const Vec3& color) { - double ambientStrength = 0.1; - double diffuseStrength = 0.8; - double specularStrength = 0.5; - int shininess = 64; - - Vec3 ambient = color * ambientStrength; - double diff = max(normal.dot(lightDir), 0.0); - Vec3 diffuse = color * (diffuseStrength * diff); - - Vec3 reflectDir = (normal * (2.0 * normal.dot(lightDir)) - lightDir).normalize(); - double spec = pow(max(viewDir.dot(reflectDir), 0.0), shininess); - Vec3 specular = Vec3(1.0, 1.0, 1.0) * (specularStrength * spec); - - return ambient + diffuse + specular; -} - -__global__ void renderKernel(unsigned char* framebuffer, Sphere* spheres, int numSpheres, Vec3 lightPos) { - int x = blockIdx.x * blockDim.x + threadIdx.x; - int y = blockIdx.y * blockDim.y + threadIdx.y; - if (x >= WIDTH || y >= HEIGHT) return; - - int pixelIndex = (y * WIDTH + x) * 3; - Vec3 rayOrigin(0, 0, 0); - Vec3 colCum(0, 0, 0); - - double spp = static_cast(SAMPLES_PER_PIXEL); - for (int sample = 0; sample < SAMPLES_PER_PIXEL; sample++) { - double u = (x + (sample / spp) - WIDTH / 2.0) / WIDTH; - double v = (y + (sample / spp) - HEIGHT / 2.0) / HEIGHT; - Vec3 rayDir(u, v, 1.0); - rayDir = rayDir.normalize(); - - for (int i = 0; i < numSpheres; ++i) { - double t; - if (spheres[i].intersect(rayOrigin, rayDir, t)) { - Vec3 hitPoint = rayOrigin + rayDir * t; - Vec3 normal = (hitPoint - spheres[i].center).normalize(); - Vec3 lightDir = (lightPos - hitPoint).normalize(); - Vec3 viewDir = -rayDir; - - colCum = colCum + phongShading(hitPoint, normal, lightDir, viewDir, spheres[i].color); - } - } + int num = 0; + for(int i = 0; i < x.size(); i++) { + if (x[i] < 1E14) std::cout << x[i] << "\n"; + if(num > 10000) break; + num++; } - // Average color across all samples - Vec3 color = colCum * (1.0 / SAMPLES_PER_PIXEL); - - framebuffer[pixelIndex] = static_cast(fmin(color.x, 1.0) * 255); - framebuffer[pixelIndex + 1] = static_cast(fmin(color.y, 1.0) * 255); - framebuffer[pixelIndex + 2] = static_cast(fmin(color.z, 1.0) * 255); -} - - - -int main() { - Sphere spheres[] = { - { Vec3(0, 0, 5), 1.0, Vec3(1.0, 0.0, 0.0) }, // Red sphere - { Vec3(-2, 1, 7), 1.0, Vec3(0.0, 1.0, 0.0) }, // Green sphere - { Vec3(2, -1, 6), 1.0, Vec3(0.0, 0.0, 1.0) } // Blue sphere - }; - int numSpheres = sizeof(spheres) / sizeof(Sphere); - Vec3 lightPos(5, 5, 0); - - unsigned char* d_framebuffer; - unsigned char* h_framebuffer = new unsigned char[WIDTH * HEIGHT * 3]; - Sphere* d_spheres; - cudaMalloc(&d_framebuffer, WIDTH * HEIGHT * 3); - cudaMalloc(&d_spheres, numSpheres * sizeof(Sphere)); - cudaMemcpy(d_spheres, spheres, numSpheres * sizeof(Sphere), cudaMemcpyHostToDevice); - - dim3 threadsPerBlock(16, 16); - dim3 numBlocks((WIDTH + threadsPerBlock.x - 1) / threadsPerBlock.x, - (HEIGHT + threadsPerBlock.y - 1) / threadsPerBlock.y); - renderKernel<<>>(d_framebuffer, d_spheres, numSpheres, lightPos); - cudaDeviceSynchronize(); - - cudaMemcpy(h_framebuffer, d_framebuffer, WIDTH * HEIGHT * 3, cudaMemcpyDeviceToHost); - saveImage("output.ppm", h_framebuffer, WIDTH, HEIGHT); - - cudaFree(d_framebuffer); - cudaFree(d_spheres); - delete[] h_framebuffer; - - std::cout << "High-resolution image saved as output.ppm" << std::endl; return 0; } diff --git a/src/objs/sphere.h b/src/objs/sphere.h deleted file mode 100644 index a9b20b4..0000000 --- a/src/objs/sphere.h +++ /dev/null @@ -1,23 +0,0 @@ -#pragma once - -#include -#include - -#include "linalg/linalg.h" - -struct Sphere { - Vec3 center; - double radius; - Vec3 color; - - __device__ bool intersect(const Vec3& rayOrigin, const Vec3& rayDir, double& t) const { - Vec3 oc = rayOrigin - center; - double b = oc.dot(rayDir); - double c = oc.dot(oc) - radius * radius; - double h = b * b - c; - if (h < 0.0) return false; - h = sqrt(h); - t = -b - h; - return true; - } -}; \ No newline at end of file diff --git a/src/tests/a.out b/src/tests/a.out deleted file mode 100755 index 3cf0e96..0000000 Binary files a/src/tests/a.out and /dev/null differ diff --git a/src/tests/deprecated_test.cu b/src/tests/deprecated_test.cu deleted file mode 100644 index b661dad..0000000 --- a/src/tests/deprecated_test.cu +++ /dev/null @@ -1,128 +0,0 @@ -#include -#include -#include -#include - -#define WIDTH 800 -#define HEIGHT 600 - -struct Vec3 { - double x, y, z; - - __host__ __device__ Vec3() : x(0), y(0), z(0) {} - __host__ __device__ Vec3(double x, double y, double z) : x(x), y(y), z(z) {} - - __host__ __device__ Vec3 operator+(const Vec3& b) const { return Vec3(x + b.x, y + b.y, z + b.z); } - __host__ __device__ Vec3 operator-(const Vec3& b) const { return Vec3(x - b.x, y - b.y, z - b.z); } - __host__ __device__ Vec3 operator*(double b) const { return Vec3(x * b, y * b, z * b); } - __host__ __device__ Vec3 operator-() const { return Vec3(-x, -y, -z); } - __host__ __device__ double dot(const Vec3& b) const { return x * b.x + y * b.y + z * b.z; } - __host__ __device__ Vec3 normalize() const { double len = sqrt(x * x + y * y + z * z); return Vec3(x / len, y / len, z / len); } -}; - -// Simple Phong lighting components -struct Sphere { - Vec3 center; - double radius; - Vec3 color; - - __device__ bool intersect(const Vec3& rayOrigin, const Vec3& rayDir, double& t) const { - Vec3 oc = rayOrigin - center; - double b = oc.dot(rayDir); - double c = oc.dot(oc) - radius * radius; - double h = b * b - c; - if (h < 0.0) return false; - h = sqrt(h); - t = -b - h; - return true; - } -}; - -__device__ Vec3 phongShading(const Vec3& point, const Vec3& normal, const Vec3& lightDir, const Vec3& viewDir, const Vec3& color) { - double ambientStrength = 0.1; - double diffuseStrength = 0.8; - double specularStrength = 0.5; - int shininess = 32; - - // Ambient - Vec3 ambient = color * ambientStrength; - - // Diffuse - double diff = max(normal.dot(lightDir), 0.0); - Vec3 diffuse = color * (diffuseStrength * diff); - - // Specular - Vec3 reflectDir = (normal * (2.0 * normal.dot(lightDir)) - lightDir).normalize(); - double spec = pow(max(viewDir.dot(reflectDir), 0.0), shininess); - Vec3 specular = Vec3(1.0, 1.0, 1.0) * (specularStrength * spec); - - return ambient + diffuse + specular; -} - -__global__ void renderKernel(unsigned char* framebuffer, Sphere sphere, Vec3 lightPos) { - int x = blockIdx.x * blockDim.x + threadIdx.x; - int y = blockIdx.y * blockDim.y + threadIdx.y; - if (x >= WIDTH || y >= HEIGHT) return; - - int pixelIndex = (y * WIDTH + x) * 3; - - Vec3 rayOrigin(0, 0, 0); - Vec3 rayDir((x - WIDTH / 2.0) / WIDTH, (y - HEIGHT / 2.0) / HEIGHT, 1.0); - rayDir = rayDir.normalize(); - - double t; - if (sphere.intersect(rayOrigin, rayDir, t)) { - Vec3 hitPoint = rayOrigin + rayDir * t; - Vec3 normal = (hitPoint - sphere.center).normalize(); - Vec3 lightDir = (lightPos - hitPoint).normalize(); - Vec3 viewDir = -rayDir; - - Vec3 color = phongShading(hitPoint, normal, lightDir, viewDir, sphere.color); - - framebuffer[pixelIndex] = static_cast(fmin(color.x, 1.0) * 255); - framebuffer[pixelIndex + 1] = static_cast(fmin(color.y, 1.0) * 255); - framebuffer[pixelIndex + 2] = static_cast(fmin(color.z, 1.0) * 255); - } else { - framebuffer[pixelIndex] = 0; - framebuffer[pixelIndex + 1] = 0; - framebuffer[pixelIndex + 2] = 0; - } -} - -void saveImage(const char* filename, unsigned char* framebuffer) { - std::ofstream imageFile(filename, std::ios::out | std::ios::binary); - imageFile << "P6\n" << WIDTH << " " << HEIGHT << "\n255\n"; - for (int i = 0; i < WIDTH * HEIGHT * 3; i++) { - imageFile << framebuffer[i]; - } - imageFile.close(); -} - -int main() { - // Initialize sphere and light source - Sphere sphere = { Vec3(0, 0, 5), 1.0, Vec3(1.0, 0.0, 0.0) }; // Red sphere - Vec3 lightPos(5, 5, 0); - - // Allocate framebuffer on device and host - unsigned char* d_framebuffer; - unsigned char* h_framebuffer = new unsigned char[WIDTH * HEIGHT * 3]; - cudaMalloc(&d_framebuffer, WIDTH * HEIGHT * 3); - - // Launch - dim3 threadsPerBlock(16, 16); - dim3 numBlocks((WIDTH + threadsPerBlock.x - 1) / threadsPerBlock.x, - (HEIGHT + threadsPerBlock.y - 1) / threadsPerBlock.y); - renderKernel<<>>(d_framebuffer, sphere, lightPos); - cudaDeviceSynchronize(); - - // Copy result back to host and save - cudaMemcpy(h_framebuffer, d_framebuffer, WIDTH * HEIGHT * 3, cudaMemcpyDeviceToHost); - saveImage("output.ppm", h_framebuffer); - - // Clean up - cudaFree(d_framebuffer); - delete[] h_framebuffer; - - std::cout << "Image saved as output.ppm" << std::endl; - return 0; -} diff --git a/src/tests/hello_world b/src/tests/hello_world deleted file mode 100755 index bb85e75..0000000 Binary files a/src/tests/hello_world and /dev/null differ diff --git a/src/tests/hello_world.cu b/src/tests/hello_world.cu deleted file mode 100644 index a3b2c6a..0000000 --- a/src/tests/hello_world.cu +++ /dev/null @@ -1,16 +0,0 @@ -#include -#include - -__global__ void hello_from_gpu() { - printf("Hello from GPU!\n"); -} - -int main() { - hello_from_gpu<<<1, 1>>>(); - - cudaDeviceSynchronize(); - - // Reset device - cudaDeviceReset(); - return 0; -} diff --git a/src/tests/output.ppm b/src/tests/output.ppm deleted file mode 100644 index 1d7b780..0000000 Binary files a/src/tests/output.ppm and /dev/null differ diff --git a/src/tests/test_heavy.cu b/src/tests/test_heavy.cu deleted file mode 100644 index 1f70108..0000000 --- a/src/tests/test_heavy.cu +++ /dev/null @@ -1,102 +0,0 @@ -#include -#include -#include -#include - -#include "linalg/linalg.h" -#include "sphere.h" -#include "img/handler.h" - -#define WIDTH 3840 -#define HEIGHT 2160 -#define SAMPLES_PER_PIXEL 8 - - -__device__ Vec3 phongShading(const Vec3& point, const Vec3& normal, const Vec3& lightDir, const Vec3& viewDir, const Vec3& color) { - double ambientStrength = 0.1; - double diffuseStrength = 0.8; - double specularStrength = 0.5; - int shininess = 64; - - Vec3 ambient = color * ambientStrength; - double diff = max(normal.dot(lightDir), 0.0); - Vec3 diffuse = color * (diffuseStrength * diff); - - Vec3 reflectDir = (normal * (2.0 * normal.dot(lightDir)) - lightDir).normalize(); - double spec = pow(max(viewDir.dot(reflectDir), 0.0), shininess); - Vec3 specular = Vec3(1.0, 1.0, 1.0) * (specularStrength * spec); - - return ambient + diffuse + specular; -} - -__global__ void renderKernel(unsigned char* framebuffer, Sphere* spheres, int numSpheres, Vec3 lightPos) { - int x = blockIdx.x * blockDim.x + threadIdx.x; - int y = blockIdx.y * blockDim.y + threadIdx.y; - if (x >= WIDTH || y >= HEIGHT) return; - - int pixelIndex = (y * WIDTH + x) * 3; - Vec3 rayOrigin(0, 0, 0); - Vec3 colCum(0, 0, 0); - - double spp = static_cast(SAMPLES_PER_PIXEL); - for (int sample = 0; sample < SAMPLES_PER_PIXEL; sample++) { - double u = (x + (sample / spp) - WIDTH / 2.0) / WIDTH; - double v = (y + (sample / spp) - HEIGHT / 2.0) / HEIGHT; - Vec3 rayDir(u, v, 1.0); - rayDir = rayDir.normalize(); - - for (int i = 0; i < numSpheres; ++i) { - double t; - if (spheres[i].intersect(rayOrigin, rayDir, t)) { - Vec3 hitPoint = rayOrigin + rayDir * t; - Vec3 normal = (hitPoint - spheres[i].center).normalize(); - Vec3 lightDir = (lightPos - hitPoint).normalize(); - Vec3 viewDir = -rayDir; - - colCum = colCum + phongShading(hitPoint, normal, lightDir, viewDir, spheres[i].color); - } - } - } - - // Average color across all samples - Vec3 color = colCum * (1.0 / SAMPLES_PER_PIXEL); - - framebuffer[pixelIndex] = static_cast(fmin(color.x, 1.0) * 255); - framebuffer[pixelIndex + 1] = static_cast(fmin(color.y, 1.0) * 255); - framebuffer[pixelIndex + 2] = static_cast(fmin(color.z, 1.0) * 255); -} - - - -int main() { - Sphere spheres[] = { - { Vec3(0, 0, 5), 1.0, Vec3(1.0, 0.0, 0.0) }, // Red sphere - { Vec3(-2, 1, 7), 1.0, Vec3(0.0, 1.0, 0.0) }, // Green sphere - { Vec3(2, -1, 6), 1.0, Vec3(0.0, 0.0, 1.0) } // Blue sphere - }; - int numSpheres = sizeof(spheres) / sizeof(Sphere); - Vec3 lightPos(5, 5, 0); - - unsigned char* d_framebuffer; - unsigned char* h_framebuffer = new unsigned char[WIDTH * HEIGHT * 3]; - Sphere* d_spheres; - cudaMalloc(&d_framebuffer, WIDTH * HEIGHT * 3); - cudaMalloc(&d_spheres, numSpheres * sizeof(Sphere)); - cudaMemcpy(d_spheres, spheres, numSpheres * sizeof(Sphere), cudaMemcpyHostToDevice); - - dim3 threadsPerBlock(16, 16); - dim3 numBlocks((WIDTH + threadsPerBlock.x - 1) / threadsPerBlock.x, - (HEIGHT + threadsPerBlock.y - 1) / threadsPerBlock.y); - renderKernel<<>>(d_framebuffer, d_spheres, numSpheres, lightPos); - cudaDeviceSynchronize(); - - cudaMemcpy(h_framebuffer, d_framebuffer, WIDTH * HEIGHT * 3, cudaMemcpyDeviceToHost); - saveImage("output.ppm", h_framebuffer, WIDTH, HEIGHT); - - cudaFree(d_framebuffer); - cudaFree(d_spheres); - delete[] h_framebuffer; - - std::cout << "High-resolution image saved as output.ppm" << std::endl; - return 0; -} diff --git a/test_read.py b/test_read.py new file mode 100644 index 0000000..bcee40c --- /dev/null +++ b/test_read.py @@ -0,0 +1,23 @@ +import numpy as np +from netCDF4 import Dataset + +# Load the NetCDF file +file_path = 'data/MERRA2_400.inst6_3d_ana_Np.20120101.nc4' +ncfile = Dataset(file_path, 'r') + +# Check the available variables in the file +print(ncfile.variables.keys()) + +U = ncfile.variables['T'][:] + +# Check the shape of the variable (it should be 3D) +print("Shape of U:", U.shape) + +# Compute the mean of the variable across all axes (for all elements in U) +U_mean = np.mean(U) + +# Print the mean +print("Mean of U:", U_mean) + +# Close the NetCDF file +ncfile.close()