Homec4science
Diffusion PHPC Shallow Water Project (master) Browse

R6678/parallel_code/compute_FV_parallel.cu
232b0a5d6bdemaster

Last Change
Raw File

compute_FV_parallel.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;
int main(){
//-------User Parameters--------------------------------------------------------
bool verbose = true; // Control messages
//-------Basic Parameters of the simulation-------------------------------------
clock_t timer = clock();
int Size = 500; // Size of map, Size*Size [km]
size_t nx = 2001; // Grid 1D size
double dx = ((double)Size)/((double)nx); // Grid spacening
float Tend = 0.2; // Simulation time in hours [hr]
int numElements = nx*nx; // Total number of elements
size_t memsize = numElements * sizeof(double); // Memory size of one array
int Ntmax = 100; // Maximum of iteration
//-------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, *d_HU, *d_HV; // Water height and x,y speeds
double *d_Ht, *d_HUt, *d_HVt; // Temporary memory of H HU and HV
double *d_Zdx, *d_Zdy; // Topology of the map
double *d_mu; // device time step
int *d_mutex;
// Tracking variables
double *dt_array; // Record the evolution time steps
string datapath = "../data/"; // Path for the data
int GPU_memory_need = 0; // Device RAM occupancy
int Host_memory_need = 0; // Hoste RAM occupancy
//-------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(Ntmax*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 device memory for computing-----------------------------------
if(verbose) cout << " Allocating device memory on host.. ";
gpuErrchk(cudaMalloc((void **) &d_H, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_HU, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_HV, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_Ht, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_HUt, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_HVt, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_Zdx, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_Zdy, memsize)); GPU_memory_need+= memsize;
gpuErrchk(cudaMalloc((void **) &d_mu, sizeof(double))); GPU_memory_need+= sizeof(double);
gpuErrchk(cudaMalloc((void **) &d_mutex, sizeof(int))); GPU_memory_need+= sizeof(int);
if(verbose) cout << GPU_memory_need/1e6 << " MB allocated\n";
//-------Set some device variables----------------------------------------------
gpuErrchk(cudaMemset(d_mu, 0.0, sizeof(double)));
gpuErrchk(cudaMemset(d_mutex, 0, sizeof(float)));
//-------Copy initial conditions from host to device----------------------------
if(verbose) cout << " Copying variables from host to device.." << endl;
gpuErrchk(cudaMemcpy(d_H, H, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HU, HU, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HV, HV, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_Zdx, Zdx, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_Zdy, Zdy, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_Ht, Ht, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HUt, HUt, memsize, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_HVt, HVt, memsize, cudaMemcpyHostToDevice));
cpy_to(Ht,H,numElements);
cpy_to(HUt,HU,numElements);
cpy_to(HVt,HV,numElements);
copy_host2device(d_H,d_HU,d_HV,H,HU,HV,memsize);
copy_host2device(d_Ht,d_HUt,d_HVt,Ht,HUt,HVt,memsize);
//-------Set grid and block dimensions------------------------------------------
int Nthreadx = 256;
int Nblockx = ((nx-2)*(nx-2) + Nthreadx -1)/Nthreadx;
dim3 threadsPerBlock(Nthreadx);
dim3 numBlocks(Nblockx);
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: " << Nblockx << endl;
cout <<"\t Nthreads \t\t\t: " << Nblockx << "x" << Nthreadx
<< " (=" << Nblockx*Nthreadx << ")" << endl;
}
//------------------------------------------------------------------------------
//-------Evolution loop------------/!\ Work in progress /!\---------------------
//------------------------------------------------------------------------------
while (T < Tend and nt < Ntmax) {
//-------Compute the time-step length-------------------------------------------
// update_dt(H,HU,HV,dt,dx,numElements);
cudaMemset(d_mu, 0.0, sizeof(double));
find_mumax_kernel<<<256,256>>>(d_H,d_HU,d_HV,d_mutex,d_mu,numElements);
cudaMemcpy(GPU_mu, d_mu, sizeof(double), cudaMemcpyDeviceToHost);
GPU_dt = dx/(sqrt(2.0)*GPU_mu[0]);
cudaDeviceSynchronize();
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------------------------------------------
cpy_to(Ht,H,numElements);
cpy_to(HUt,HU,numElements);
cpy_to(HVt,HV,numElements);
// copy_temp_variables_kernel<<<numBlocks,threadsPerBlock>>>(d_H,d_HU,d_HV,d_Ht,d_HUt,d_HVt,nx);
// cudaDeviceSynchronize();
//-------Enforce boundary conditions--------------------------------------------
enforce_BC(Ht, HUt, HVt, nx);
//dim3 numblock(16,4); //(four times 16*128 = 2048 approx nx)
// enforce_BC_kernel<<<numblock,128>>>(d_Ht,d_HUt,d_HVt,nx);
// cudaDeviceSynchronize();
//-------Compute a time-step----------------------------------------------------
C = (.5*GPU_dt/dx);
//FV_time_step(H,HU,HV,Zdx,Zdy,Ht,HUt,HVt,C,dt[0],nx);
//FV_time_step_kernel<<<numBlocks,threadsPerBlock>>>(d_H,d_HU,d_HV,d_Zdx,d_Zdy,d_Ht,d_HUt,d_HVt,C,dt[0],nx);
copy_host2device(d_Ht,d_HUt,d_HVt,Ht,HUt,HVt,memsize);
cudaDeviceSynchronize();
FV_iterator_kernel<<<numBlocks,threadsPerBlock>>>(d_H,d_HU,d_HV,d_Zdx,d_Zdy,d_Ht,d_HUt,d_HVt,C,GPU_dt,nx);
cudaDeviceSynchronize();
copy_device2host(H,HU,HV,d_H,d_HU,d_HV,memsize);
//-------Impose tolerances------------------------------------------------------
impose_tolerances(H,HU,HV,numElements);
//impose_tolerances_kernel<<<numBlocks,threadsPerBlock>>>(d_H,d_HU,d_HV,numElements);
//cudaDeviceSynchronize();
//-------Store time step evolution and update time and counter------------------
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;
}
cudaMemcpy(H, d_H, memsize, cudaMemcpyDeviceToHost);
//-------Save solution to disk--------------------------------------------------
ostringstream soutfilename;
soutfilename <<"../output/CUDA_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/CUDA_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]), Ntmax*sizeof(double));
fout.close();
//-------Free device global memory----------------------------------------------
if(verbose) cout << " Free device memory space.." << endl;
cudaFree(d_H); cudaFree(d_HU); cudaFree(d_HV); cudaFree(d_Zdx);
cudaFree(d_Zdy); cudaFree(d_Ht); cudaFree(d_HUt); cudaFree(d_HVt);
cudaFree(d_mu); cudaFree(d_mutex);
//-------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;
timer = (double)(timer)/CLOCKS_PER_SEC*1000;
if(verbose){
cout << "Ellapsed time : " << timer/60000 << "min "
<< (timer/1000)%60 << "s " << timer%1000 << "ms" << endl;
}
return 0;
}
c4science · Help