Homec4science
Diffusion PHPC Shallow Water Project (master) Browse

R6678/multi_gpu_code/compute_parallel_multi_GPU.cu
b053dad92701master

Last Change
Raw File

compute_parallel_multi_GPU.cu
View Options

/*
Parallel version in CUDA/C++ of compute.cpp code
*/
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <cstring>
#include <sstream>
#include "kernels.cuh" //parallel kernels
#include "functions.h" //Sequential functions
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true)
{
if (code != cudaSuccess)
{
fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
if (abort) getchar();
}
}
using namespace std;
void print_device_properties(int dev_id){
cudaDeviceProp prop;
gpuErrchk(cudaGetDeviceProperties(&prop,dev_id));
cout << prop.name << " (" << prop.totalGlobalMem/1e9 <<"GB RAM, id = "<<dev_id<<")";
}
int main(int argc, char *argv[]){
//-------User Parameters--------------------------------------------------------
bool verbose = true; // Control messages
int Nthreads = 256; // Number of threads per block
int nx = 2001; // Grid 1D size (2001,4001,8001)
int Ntmax = 100; // Maximum of iteration, =0 for only Tend bound
float Tend = 0.2; // Simulation time in hours [hr]
string author = "CUDA";
int num_devices; // Number of GPUs in use
gpuErrchk(cudaGetDeviceCount(&num_devices)); // Default is the detected number of device
if(argc>1){
string filename = argv[1];
int *NX = (int *)malloc(sizeof(int));
int *NTHREADS = (int *)malloc(sizeof(int));
int *VERBOSE = (int *)malloc(sizeof(int));
int *NTMAX = (int *)malloc(sizeof(int));
float *TEND = (float *)malloc(sizeof(float));
int *AUTHOR = (int *)malloc(sizeof(int));
int * NDEVICES = (int *)malloc(sizeof(int));
read_parameters(filename,NX,NTHREADS,VERBOSE,NTMAX,TEND,AUTHOR,NDEVICES);
verbose = *VERBOSE;
Nthreads = *NTHREADS;
nx = *NX;
Ntmax = *NTMAX;
Tend = *TEND;
if(*AUTHOR == 0) author = "GTX";
if(*AUTHOR == 1) author = "Deneb";
if(num_devices >= *NDEVICES){
if(verbose) cout<<"#Detected Devices : "<< num_devices << endl;
if(verbose) cout<<"#Set Devices : "<< *NDEVICES << endl;
num_devices = *NDEVICES;
}
else{
cerr<<"ERROR : #set devices("<<NDEVICES<<") > #detected devices ("<<num_devices<<")\n";
return 0;
}
free(NX); free(NTHREADS); free(VERBOSE); free(NTMAX); free(AUTHOR); free(NDEVICES);
}
//-------Basic Parameters of the simulation-------------------------------------
clock_t timer = clock(); // Timer for performance measurement
unsigned long int ops; // Total number of operation
double gflops = 0; // Performance in Gflops
int Size = 500; // Size of map, Size*Size [km]
double dx = ((double)Size)/((double)nx); // Grid spacening
int numElements = nx*nx; // Total number of elements
size_t memsize = numElements * sizeof(double); // Memory size of one array
//-------Simulation variables HOST----------------------------------------------
double T = 0.0; // Time
int nt = 0; // Iteration counter
double C = 0.0; // Coefficient 1/2*dt/dx
double *H, *HU, *HV; // Water height and x,y speeds
double *Ht, *HUt, *HVt; // Temporary memory of H HU and HV
double *Zdx, *Zdy; // Topology of the map
double *dt; // Host Time step
double *GPU_mu;
double GPU_dt;
//-------Simulation variables DEVICE--------------------------------------------
double *d_H[num_devices]; // Water height
double *d_HU[num_devices], *d_HV[num_devices]; // x,y Water velocities
double *d_Ht[num_devices]; // Temporary memory of H
double *d_HTt[num_devices], *d_HVt[num_devices]; // Temporary memory of HU and HV
double *d_Zdx[num_devices], *d_Zdy[num_devices]; // Topology of the map
double *d_mu[num_devices]; // device time step
int *d_mutex[num_devices];
// Tracking variables
double *dt_array; // Record the evolution time steps
string datapath = "../data/"; // Path for the data
unsigned long int GPU_memory_need = 0; // Device RAM occupancy
unsigned long int Host_memory_need = 0; // Hoste RAM occupancy
//-------MULTI GPU Variables----------------------------------------------------
int startMU; int endMU;
int startCPY; int endCPY;
int startBC; int endBC;
int startFV; int endFV;
int startTOL; int endTOL;
int chunk_size = nx/num_devices * nx;
double chunk_memsize = chunk_size * sizeof(double);
//-------Allocate host memory for loading the initial conditions----------------
if(verbose) cout << " Allocating host memory .. ";
H = (double *)malloc(memsize); Host_memory_need += memsize;
HU = (double *)malloc(memsize); Host_memory_need += memsize;
HV = (double *)malloc(memsize); Host_memory_need += memsize;
Ht = (double *)malloc(memsize); Host_memory_need += memsize;
HUt = (double *)malloc(memsize); Host_memory_need += memsize;
HVt = (double *)malloc(memsize); Host_memory_need += memsize;
Zdx = (double *)malloc(memsize); Host_memory_need += memsize;
Zdy = (double *)malloc(memsize); Host_memory_need += memsize;
dt = (double *)malloc(sizeof(double)); Host_memory_need += sizeof(double);
GPU_mu = (double *)malloc(sizeof(double)); Host_memory_need += sizeof(double);
dt_array = (double *)malloc(fmaxf(Ntmax,100)*sizeof(double)); Host_memory_need += Ntmax*sizeof(double);
if(verbose) cout << Host_memory_need/1e6 << " MB allocated\n";
//-------Load initial state on host memory--------------------------------------
load_initial_data(H,HU,HV,Zdx,Zdy,datapath,nx,Size,Tend,numElements,verbose);
//-------Allocate devices memory for computing----------------------------------
if(verbose) cout << " Allocating devices memory.."<<endl;
for(int dev_id = 0; dev_id < num_devices; ++dev_id){
chunk_memsize = chunk_size * sizeof(double);
if(dev_id = num_devices-1) chunk_memsize += nx%num_devices * nx * sizeof(double);
gpuErrchk(cudaSetDevice( dev_id ));
gpuErrchk(cudaFree(0));
gpuErrchk(cudaMalloc((void **) d_H+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_HU+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_HV+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_Ht+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_HTt+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_HVt+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_Zdx+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_Zdy+dev_id, chunk_memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) d_mu+dev_id, sizeof(double))); GPU_memory_need+= sizeof(double);
gpuErrchk(cudaMalloc((void **) d_mutex+dev_id, sizeof(int))); GPU_memory_need+= sizeof(int);
if(verbose) cout << GPU_memory_need/1e6 << " MB allocated on ";print_device_properties(dev_id);
//-------Set some device variables----------------------------------------------
gpuErrchk(cudaMemset(d_mu[dev_id], 0.0, sizeof(double)));
gpuErrchk(cudaMemset(d_mutex[dev_id], 0, sizeof(float)));
//-------Copy initial conditions from host to device----------------------------
if(verbose) cout << " Copying variables from host to device " << dev_id << endl;
gpuErrchk(cudaMemcpy(d_H[dev_id], H +dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HU[dev_id], HU +dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HV[dev_id], HV +dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_Zdx[dev_id], Zdx+dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_Zdy[dev_id], Zdy+dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_Ht[dev_id], Ht +dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HTt[dev_id], HUt+dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HVt[dev_id], HVt+dev_id*chunk_size, chunk_memsize, cudaMemcpyHostToDevice));
}
//-------Set grid and block dimensions------------------------------------------
// For finite volume kernel
int NblocksFV = ((nx-2)*(nx-2) + Nthreads -1)/Nthreads;
// For tolerances kernel
int NblocksBC = pow(2,ceil(log2(nx)))/Nthreads;
dim3 GridDimBC(NblocksBC,4);
// For Enforce BC kernel
int NblocksTol = (nx*nx + Nthreads -1)/Nthreads;
if(verbose){
cout << "parallel FV iterator kernel :" << endl;
cout <<"\t Number of elements \t\t: " << (nx-2)*(nx-2) << endl;
cout <<"\t Number of blocks needed \t: " << NblocksFV << endl;
cout <<"\t Nthreads \t\t\t: " << NblocksFV << "x" << Nthreads
<< " (=" << NblocksFV*Nthreads << ")" << endl;
cout << "parallel Enforce BC kernel :" << endl;
cout <<"\t Number of elements \t\t: " << 4*nx << endl;
cout <<"\t Number of blocks needed \t: " << NblocksBC << endl;
cout <<"\t Nthreads \t\t\t: " << NblocksBC << "x4x" << Nthreads
<< " (=" << NblocksBC*Nthreads*4 << ")" << endl;
cout << "parallel Tolerances kernel :" << endl;
cout <<"\t Number of elements \t\t: " << nx*nx << endl;
cout <<"\t Number of blocks needed \t: " << NblocksTol << endl;
cout <<"\t Nthreads \t\t\t: " << NblocksTol<< "x" << Nthreads
<< " (=" << NblocksTol*Nthreads << ")" << endl;
}
//------------------------------------------------------------------------------
//-------Evolution loop------------/!\ Work in progress /!\---------------------
//------------------------------------------------------------------------------
while (T < Tend and (nt < Ntmax or Ntmax == 0)) {
//-------Compute the time-step length-------------------------------------------
for ( int dev_id = 0; dev_id < num_devices; ++dev_id){
startMU = dev_id * chunk_size;
endMU = (dev_id+1)* chunk_size;
gpuErrchk(cudaSetDevice(dev_id));
find_mumax_kernel<<<Nthreads,Nthreads>>>
(d_H[dev_id],d_HU[dev_id],d_HV[dev_id],d_mutex[dev_id],d_mu[dev_id],startMU,endMU,nx);
gpuErrchk(cudaMemcpy(GPU_mu, d_mu[dev_id], sizeof(double), cudaMemcpyDeviceToHost));
}
GPU_dt = dx/(sqrt(2.0)*GPU_mu[0]);
if(T+GPU_dt > Tend){
GPU_dt = Tend-T;
}
//-------Print status-----------------------------------------------------------
if(verbose){
cout << " Computing for T=" << T+GPU_dt
<< "\t("<< 100*(T+GPU_dt)/Tend << "%),\t"
<< "dt=" << GPU_dt << "\n";
// cout << " CPU dt \t : " << dt[0] << "\n";
cout << " GPU dt \t : " << GPU_dt << "\n";
}
//-------Copy solution to temp storage------------------------------------------
for(int dev_id = 0; dev_id < num_devices; ++dev_id){
startCPY = dev_id * chunk_size;
endCPY = (dev_id+1)* chunk_size;
gpuErrchk(cudaSetDevice(dev_id));
copy_temp_variables_kernel<<<NblocksTol,Nthreads>>>
(d_H[dev_id],d_HU[dev_id],d_HV[dev_id],d_Ht[dev_id],
d_HTt[dev_id],d_HVt[dev_id],startCPY,endCPY,nx);
}
cudaDeviceSynchronize();
//-------Enforce boundary conditions--------------------------------------------
for(int dev_id = 0; dev_id < num_devices; ++dev_id){
startBC = dev_id * chunk_size;
endBC = (dev_id+1)* chunk_size;
gpuErrchk(cudaSetDevice(dev_id));
enforce_BC_kernel<<<GridDimBC,Nthreads>>>
(d_Ht[dev_id],d_HTt[dev_id],d_HVt[dev_id],startBC,endBC,nx);
}
cudaDeviceSynchronize();
//-------Compute a time-step----------------------------------------------------
C = (.5*GPU_dt/dx);
for(int dev_id = 0; dev_id < num_devices; ++dev_id){
startFV = dev_id * chunk_size;
endFV = (dev_id+1)* chunk_size;
gpuErrchk(cudaSetDevice(dev_id));
FV_iterator_kernel<<<NblocksFV,Nthreads>>>
(d_H[dev_id],d_HU[dev_id],d_HV[dev_id],d_Zdx[dev_id],d_Zdy[dev_id],
d_Ht[dev_id],d_HTt[dev_id],d_HVt[dev_id],C,GPU_dt,startFV,endFV,nx);
}
cudaDeviceSynchronize();
//-------Impose tolerances------------------------------------------------------
for(int dev_id = 0; dev_id < num_devices; ++dev_id){
startTOL = dev_id * chunk_size;
endTOL = (dev_id+1)* chunk_size;
gpuErrchk(cudaSetDevice(dev_id));
impose_tolerances_kernel<<<NblocksTol,Nthreads>>>
(d_H[dev_id],d_HU[dev_id],d_HV[dev_id],startTOL,endTOL,nx);
}
//-------Store time step evolution and update time and counter------------------
if(nt < fmaxf(100,Ntmax)) dt_array[nt] = GPU_dt;
T = T + GPU_dt;
nt++;
cudaDeviceSynchronize();
}
//------------------------------------------------------------------------------
//-------Copy device result to the host memory----------------------------------
if(verbose){
cout << " Copy the output data from the CUDA device to the host memory" << endl;
}
gpuErrchk(cudaMemcpy(H, d_H[0], memsize, cudaMemcpyDeviceToHost));
//-------Save solution to disk--------------------------------------------------
ostringstream soutfilename;
soutfilename <<"../output/"<<author<<"_Solution_nx"<<to_string(nx)<<"_"<<to_string(Size)<<"km_T"<<Tend<<"_h.bin"<< setprecision(2);
string outfilename = soutfilename.str();
ofstream fout;
fout.open(outfilename, std::ios::out | std::ios::binary);
if(verbose) cout<<" Writing solution in "<<outfilename<<endl;
fout.write(reinterpret_cast<char*>(&H[0]), numElements*sizeof(double));
fout.close();
//-------Save time step history-------------------------------------------------
ostringstream soutfilename2;
soutfilename2 <<"../output/"<<author<<"_dt_nx"<<to_string(nx)<<"_"<<to_string(Size)<<"km_T"<<Tend<<"_h.bin"<< setprecision(2);
outfilename = soutfilename2.str();
fout.open(outfilename, std::ios::out | std::ios::binary);
if(verbose) cout<<" Writing solution in "<<outfilename<<endl;
fout.write(reinterpret_cast<char*>(&dt_array[0]), fmaxf(100,Ntmax)*sizeof(double));
fout.close();
//-------Free device global memory----------------------------------------------
if(verbose) cout << " Free device memory space.." << endl;
cudaFree(d_H[0]); cudaFree(d_HU[0]); cudaFree(d_HV[0]); cudaFree(d_Zdx[0]);
cudaFree(d_Zdy[0]); cudaFree(d_Ht[0]); cudaFree(d_HTt[0]); cudaFree(d_HVt[0]);
cudaFree(d_mu[0]); cudaFree(d_mutex[0]);
//-------Free host memory-------------------------------------------------------
if(verbose) cout << " Free host memory space.." << endl;
free(H); free(HU); free(HV); free(Zdx); free(Zdy);
free(Ht); free(HUt); free(HVt); free(dt_array); free(dt); free(GPU_mu);
//-------End of the timer-----------------------------------------------------
timer = clock()-timer; // in ticks
ops = nt*( 15 + 2 + 11 + 30 + 30 + 1 )*pow(nx,2);
gflops = ops/1e9/(timer/CLOCKS_PER_SEC);
timer = (double)(timer)/CLOCKS_PER_SEC*1000; // ticks to ms
if(verbose){
cout << "Ellapsed time : " << timer/60000 << "min "
<< (timer/1000)%60 << "s " << timer%1000 << "ms" << endl;
cout << "Average performance : " << gflops << "gflops" << endl;
}
else{cout << nx << "," << Nthreads << "," << Ntmax << "," << timer/1000.0 << endl;}
return 0;
}
c4science · Help