From 0db764c52f3cef16a9e1ea0ff5a4422a565f4d81 Mon Sep 17 00:00:00 2001 From: Robin Date: Fri, 20 Dec 2024 18:46:18 +0100 Subject: [PATCH] FieldData structs implemented --- load-modules.sh | 2 +- src/hurricanedata/datareader.cu | 60 ------------------------- src/hurricanedata/datareader.h | 10 ----- src/hurricanedata/fielddata.cu | 12 +++++ src/hurricanedata/fielddata.h | 44 ++++++++++++++++++ src/hurricanedata/gpubuffer.cu | 70 +++++++++++++++++++++++++++++ src/hurricanedata/gpubuffer.h | 24 ++++++++++ src/hurricanedata/gpubufferhelper.h | 10 +++++ src/main.cu | 33 +++++++++----- test_read.py | 21 +++++---- 10 files changed, 196 insertions(+), 90 deletions(-) delete mode 100644 src/hurricanedata/datareader.cu delete mode 100644 src/hurricanedata/datareader.h create mode 100644 src/hurricanedata/fielddata.cu create mode 100644 src/hurricanedata/fielddata.h create mode 100644 src/hurricanedata/gpubuffer.cu create mode 100644 src/hurricanedata/gpubuffer.h create mode 100644 src/hurricanedata/gpubufferhelper.h diff --git a/load-modules.sh b/load-modules.sh index da7a5c5..c6fcab4 100755 --- a/load-modules.sh +++ b/load-modules.sh @@ -3,4 +3,4 @@ # Load the needed modules in Habrok ml CUDA -ml netCDF-C++4 \ No newline at end of file +ml netCDF-C++4 diff --git a/src/hurricanedata/datareader.cu b/src/hurricanedata/datareader.cu deleted file mode 100644 index 62dc72b..0000000 --- a/src/hurricanedata/datareader.cu +++ /dev/null @@ -1,60 +0,0 @@ -#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; -} - -std::pair loadDataToDevice(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(); - } - - // Store NetCDF variable in pinned memory on host - float *h_array; - - cudaMallocHost(&h_array, sizeof(float)*length); - - var.getVar(h_array); - - // Copy data to device - float *d_array; - - cudaError_t status = cudaMalloc(&d_array, sizeof(float)*length); - if (status != cudaSuccess) - cout << "Error allocating memory: " << status << "\n"; - - cudaMemcpyAsync(d_array, h_array, sizeof(float)*length, cudaMemcpyHostToDevice); - - cudaDeviceSynchronize(); // Heavy hammer synchronisation // TODO: Use streams - - cudaFreeHost(h_array); - - return std::pair(d_array, length); -} \ No newline at end of file diff --git a/src/hurricanedata/datareader.h b/src/hurricanedata/datareader.h deleted file mode 100644 index 75a7e53..0000000 --- a/src/hurricanedata/datareader.h +++ /dev/null @@ -1,10 +0,0 @@ -#ifndef DATAREADER_H -#define DATAREADER_H - -#include -#include - -std::vector readData(std::string path, std::string variableName); -std::pair loadDataToDevice(std::string path, std::string variableName); - -#endif //DATAREADER_H diff --git a/src/hurricanedata/fielddata.cu b/src/hurricanedata/fielddata.cu new file mode 100644 index 0000000..3683275 --- /dev/null +++ b/src/hurricanedata/fielddata.cu @@ -0,0 +1,12 @@ +#include "fielddata.h" + +// __device__ float getVal( +// const FieldMetadata &md, +// const FieldData &d, +// const size_t &timeInd, +// const size_t &lonInd, +// const size_t &latInd, +// const size_t &levInd +// ) { +// return d.valArrays[0][timeInd]; +// } \ No newline at end of file diff --git a/src/hurricanedata/fielddata.h b/src/hurricanedata/fielddata.h new file mode 100644 index 0000000..fac8713 --- /dev/null +++ b/src/hurricanedata/fielddata.h @@ -0,0 +1,44 @@ +#ifndef FIELDDATA_H +#define FIELDDATA_H + +#include + +struct FieldMetadata { + size_t widthSize; // Number of different longitudes + size_t heightSize; // Number of different latitudes + size_t depthSize; // Number of different levels + + // lons is a managed Unified Memory array of size widthCount that indicates + // that getVal(t, i, j, k) is a value with longitude of lons[i]. + // The other such arrays are similarly defined. + double *lons; + double *lats; + double *levs; +}; + +using FieldMetadata = FieldMetadata; + +struct FieldData { + static constexpr size_t FILESNUM = 2; // Number of files stored in a FieldData struct. + + // An array of length FILESNUM storing pointers to 4D arrays stored in device memory. + float *valArrays[FILESNUM]; + + size_t timeSize; // Number of different times + // times is a managed Unified Memory array of size timeSize that indicates + // that getVal(md, d, t, i, j, k) is a value at time times[t]. + int *times; +}; + +using FieldData = FieldData; + +// __device__ float getVal( +// const FieldMetadata &md, +// const FieldData &d, +// const size_t &timeInd, +// const size_t &lonInd, +// const size_t &latInd, +// const size_t &levInd +// ); + +#endif //FIELDDATA_H diff --git a/src/hurricanedata/gpubuffer.cu b/src/hurricanedata/gpubuffer.cu new file mode 100644 index 0000000..d0a8946 --- /dev/null +++ b/src/hurricanedata/gpubuffer.cu @@ -0,0 +1,70 @@ +#include "gpubuffer.h" +#include "fielddata.h" +#include "gpubufferhelper.h" + +#include + + +using namespace std; +using namespace netCDF; + +GPUBuffer::GPUBuffer(std::string path, std::string variableName): +path(path), variableName(variableName) { + NcFile data(path, NcFile::read); + + multimap vars = data.getVars(); + + cudaMallocManaged(&fmd, sizeof(FieldMetadata)); + + readAndAllocateAxis(&fmd->lons, &fmd->widthSize, vars.find("lon")->second); + readAndAllocateAxis(&fmd->lats, &fmd->heightSize, vars.find("lat")->second); + readAndAllocateAxis(&fmd->levs, &fmd->depthSize, vars.find("lev")->second); +} + + + +FieldData GPUBuffer::nextFieldData() { + NcFile data(path, NcFile::read); + + multimap vars = data.getVars(); + + FieldData fd; + size_t timeSize; + readAndAllocateAxis(&fd.times, &fd.timeSize, vars.find("time")->second); + + NcVar var = vars.find(variableName)->second; + + int length = 1; + for (NcDim dim: var.getDims()) { + length *= dim.getSize(); + } + + // Store NetCDF variable in pinned memory on host + float *h_array; + + cudaMallocHost(&h_array, sizeof(float)*length); + + var.getVar(h_array); + + // Copy data to device + // float *d_array; + + cudaError_t status = cudaMalloc(&fd.valArrays[0], sizeof(float)*length); + if (status != cudaSuccess) + cout << "Error allocating memory: " << status << "\n"; + + cudaMemcpyAsync(fd.valArrays[0], h_array, sizeof(float)*length, cudaMemcpyHostToDevice); + + cudaDeviceSynchronize(); // Heavy hammer synchronisation // TODO: Use streams + + cudaFreeHost(h_array); + + return fd; +} + +GPUBuffer::~GPUBuffer() { + cudaFree(fmd->lons); + cudaFree(fmd->lats); + cudaFree(fmd->levs); + cudaFree(fmd); +} \ No newline at end of file diff --git a/src/hurricanedata/gpubuffer.h b/src/hurricanedata/gpubuffer.h new file mode 100644 index 0000000..06ee9c6 --- /dev/null +++ b/src/hurricanedata/gpubuffer.h @@ -0,0 +1,24 @@ +#ifndef GPUBUFFER_H +#define GPUBUFFER_H + +#include "fielddata.h" + +#include + +class GPUBuffer { +public: + GPUBuffer(std::string path, std::string variableName); + + FieldData nextFieldData(); + + ~GPUBuffer(); + + FieldMetadata *fmd; +private: + // TODO: Implement GPUBuffer + std::string path; + std::string variableName; + +}; + +#endif //GPUBUFFER_H diff --git a/src/hurricanedata/gpubufferhelper.h b/src/hurricanedata/gpubufferhelper.h new file mode 100644 index 0000000..2f54634 --- /dev/null +++ b/src/hurricanedata/gpubufferhelper.h @@ -0,0 +1,10 @@ +#include +#include +template +void readAndAllocateAxis(T** axis_ptr, size_t *size, const netCDF::NcVar &var) { + assert(var.getDimCount() == 1); + netCDF::NcDim dim = var.getDim(0); + *size = dim.getSize(); + cudaError_t status = cudaMallocManaged(axis_ptr, *size*sizeof(T)); + var.getVar(*axis_ptr); +} \ No newline at end of file diff --git a/src/main.cu b/src/main.cu index 1f5e837..9117e63 100644 --- a/src/main.cu +++ b/src/main.cu @@ -1,19 +1,31 @@ -#include "hurricanedata/datareader.h" +#include "hurricanedata/gpubuffer.h" #include #include #include #include +__device__ float getVal( + const FieldMetadata &md, + const FieldData &d, + const size_t &timeInd, + const size_t &lonInd, + const size_t &latInd, + const size_t &levInd +) { + return d.valArrays[0][timeInd]; +} + // Not parallel computation -__global__ void computeMean(float *ans, size_t *masked_vals, size_t n, float *x) { +__global__ void computeMean(float *ans, size_t *masked_vals, const FieldMetadata &fmd, FieldData fd) { float sum = 0; size_t num_not_masked_values = 0; size_t num_masked_values = 0; - for (int i = 0; i < n; i++) { - if (x[i] < 1E14) { /* If x is not missing value */ + for (int i = 0; i < fmd.widthSize*fmd.heightSize*fmd.depthSize*fd.timeSize; i++) { + double xi = getVal(fmd, fd, i, 0, 0, 0); + if (xi < 1E14) { /* If x is not missing value */ num_not_masked_values++; - sum += x[i]; + sum += xi; } else { num_masked_values++; } @@ -25,7 +37,9 @@ __global__ void computeMean(float *ans, size_t *masked_vals, size_t n, float *x) int main() { std::string path = "data/MERRA2_400.inst6_3d_ana_Np.20120101.nc4"; std::string variable = "T"; - auto arr = loadDataToDevice(path, variable); + GPUBuffer buffer{path, variable}; + + auto fd = buffer.nextFieldData(); float *ptr_mean; cudaMallocManaged(&ptr_mean, sizeof(float)); @@ -33,15 +47,14 @@ int main() { size_t *ptr_masked; cudaMallocManaged(&ptr_masked, sizeof(size_t)); - computeMean<<<1, 1>>>(ptr_mean, ptr_masked, arr.second, arr.first); + computeMean<<<1, 1>>>(ptr_mean, ptr_masked, *buffer.fmd, fd); cudaDeviceSynchronize(); - std::cout << "Mean = " << *ptr_mean << " calculated from " << arr.second << " values where " << *ptr_masked << " are masked values.\n"; + std::cout << "Mean = " << *ptr_mean << " values where " << *ptr_masked << " are masked values.\n"; - cudaFree(arr.first); + cudaFree(fd.valArrays[0]); cudaFree(ptr_mean); cudaFree(ptr_masked); - return 0; } diff --git a/test_read.py b/test_read.py index 90eaa69..19bbcfc 100644 --- a/test_read.py +++ b/test_read.py @@ -1,6 +1,5 @@ import numpy as np from netCDF4 import Dataset -from math import prod # Load the NetCDF file file_path = 'data/MERRA2_400.inst6_3d_ana_Np.20120101.nc4' @@ -11,20 +10,24 @@ print(ncfile.variables.keys()) U = ncfile.variables['T'][:] -# Check the shape of the variable -print(f"Shape of U: {U.shape} and total length is {prod(U.shape)}") +# Check the shape of the variable +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) +U_sum = np.sum(U) # Print the mean print("Mean of U:", U_mean) -print(f"{U[0,0,0,1]=}") -is_masked = np.ma.isMaskedArray(U) -print(f"Is U a masked array? {is_masked}") -masked_count = np.ma.count_masked(U) -print("Number of masked values in U:", masked_count) +print("Sum of U:", U_sum) + +sumval = 0 +row = U[0,0,100] +for val in row: + if not np.ma.is_masked(val): + sumval += val +print(f"Why does {np.sum(row)=} not equal {sumval=} ?!") # Close the NetCDF file -ncfile.close() +ncfile.close() \ No newline at end of file