loaded data using netcdfcxx

.gitignore

@@ -0,0 +1,3 @@
+build
+.vscode
+data

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

src/hurricanedata/datareader.cu

@@ -0,0 +1,25 @@
+#include "datareader.h"
+
+#include <netcdf>
+
+using namespace std;
+using namespace netCDF;
+
+std::vector<float> readData(std::string path, std::string variableName) {
+    netCDF::NcFile data(path, netCDF::NcFile::read);
+
+    multimap<string, NcVar> vars = data.getVars();
+
+    NcVar var = vars.find(variableName)->second;   
+
+    int length = 1;
+    for (NcDim dim: var.getDims()) {
+        length *= dim.getSize();
+    }
+
+    vector<float> vec(length);
+
+    var.getVar(vec.data());
+
+    return vec;
+}

src/hurricanedata/datareader.h

@@ -0,0 +1,9 @@
+#ifndef DATAREADER_H
+#define DATAREADER_H
+
+#include <vector>
+#include <string>
+
+std::vector<float> readData(std::string path, std::string variableName);
+
+#endif //DATAREADER_H

src/img/handler.h

@@ -1,13 +0,0 @@
-#pragma once
-
-#include <fstream>
-
-
-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();
-}

src/linalg/linalg.h

@@ -1,9 +0,0 @@
-#pragma once
-
-#ifndef LINALG_H
-#define LINALG_H
-
-#include "vec.h"
-#include "mat.h"
-
-#endif // LINALG_H

src/linalg/mat.h

@@ -1 +0,0 @@
-#pragma once

src/linalg/vec.h

@@ -1,18 +0,0 @@
-#pragma once
-
-#include <cuda_runtime.h>
-#include <cmath>
-
-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); }
-};

src/main.cu

@@ -1,102 +1,21 @@
+#include "hurricanedata/datareader.h"
+
 #include <cuda_runtime.h>
 #include <device_launch_parameters.h>
 #include <iostream>
 #include <cmath>
 
-#include "linalg/linalg.h"
-#include "objs/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<double>(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<unsigned char>(fmin(color.x, 1.0) * 255);
-    framebuffer[pixelIndex + 1] = static_cast<unsigned char>(fmin(color.y, 1.0) * 255);
-    framebuffer[pixelIndex + 2] = static_cast<unsigned char>(fmin(color.z, 1.0) * 255);
-}
-
-
-
 int main() {
-    Sphere spheres[] = {
+    std::string path = "data/MERRA2_400.inst6_3d_ana_Np.20120101.nc4";
-        { Vec3(0, 0, 5), 1.0, Vec3(1.0, 0.0, 0.0) },  // Red sphere
+    std::string variable = "U";
-        { Vec3(-2, 1, 7), 1.0, Vec3(0.0, 1.0, 0.0) }, // Green sphere
+    auto x = readData(path, variable);
-        { 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;
+    int num = 0;
-    unsigned char* h_framebuffer = new unsigned char[WIDTH * HEIGHT * 3];
+    for(int i = 0; i < x.size(); i++) {
-    Sphere* d_spheres;
+        if (x[i] < 1E14) std::cout << x[i] << "\n";
-    cudaMalloc(&d_framebuffer, WIDTH * HEIGHT * 3);
+        if(num > 10000) break;
-    cudaMalloc(&d_spheres, numSpheres * sizeof(Sphere));
+        num++;
-    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<<<numBlocks, threadsPerBlock>>>(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;
 }

src/objs/sphere.h

@@ -1,23 +0,0 @@
-#pragma once
-
-#include <cuda_runtime.h>
-#include <cmath>
-
-#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;
-    }
-};

src/tests/deprecated_test.cu

@@ -1,128 +0,0 @@
-#include <cuda_runtime.h>
-#include <device_launch_parameters.h>
-#include <iostream>
-#include <fstream>
-
-#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<unsigned char>(fmin(color.x, 1.0) * 255);
-        framebuffer[pixelIndex + 1] = static_cast<unsigned char>(fmin(color.y, 1.0) * 255);
-        framebuffer[pixelIndex + 2] = static_cast<unsigned char>(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<<<numBlocks, threadsPerBlock>>>(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;
-}

src/tests/hello_world.cu

@@ -1,16 +0,0 @@
-#include <iostream>
-#include <cuda_runtime.h>
-
-__global__ void hello_from_gpu() {
-    printf("Hello from GPU!\n");
-}
-
-int main() {
-    hello_from_gpu<<<1, 1>>>();
-
-    cudaDeviceSynchronize();
-
-    // Reset device
-    cudaDeviceReset();
-    return 0;
-}

src/tests/test_heavy.cu

@@ -1,102 +0,0 @@
-#include <cuda_runtime.h>
-#include <device_launch_parameters.h>
-#include <iostream>
-#include <cmath>
-
-#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<double>(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<unsigned char>(fmin(color.x, 1.0) * 255);
-    framebuffer[pixelIndex + 1] = static_cast<unsigned char>(fmin(color.y, 1.0) * 255);
-    framebuffer[pixelIndex + 2] = static_cast<unsigned char>(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<<<numBlocks, threadsPerBlock>>>(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;
-}

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()