diff --git a/src/force.cu b/src/force.cu index a04f3bb..d68567b 100644 --- a/src/force.cu +++ b/src/force.cu @@ -99,7 +99,9 @@ __global__ void calc_force( extern "C" { bool initialized = false; -Atom c_atom; +static Atom c_atom; +int *c_neighs; +int *c_neigh_numneigh; double computeForce( bool reneighbourHappenend, @@ -117,9 +119,12 @@ double computeForce( MD_FLOAT sigma6 = param->sigma6; MD_FLOAT epsilon = param->epsilon; #endif + printf("-1\r\n"); cudaProfilerStart(); + printf("0\r\n"); + for(int i = 0; i < Nlocal; i++) { fx[i] = 0.0; fy[i] = 0.0; @@ -140,6 +145,8 @@ double computeForce( c_atom.Nmax = atom->Nmax; c_atom.ntypes = atom->ntypes; + printf("0.1\r\n"); + /* int nDevices; cudaGetDeviceCount(&nDevices); @@ -160,6 +167,8 @@ double computeForce( // HINT: Run with cuda-memcheck ./MDBench-NVCC in case of error // HINT: Only works for data layout = AOS!!! + printf("1\r\n"); + if(!initialized) { checkCUDAError( "c_atom.x malloc", cudaMalloc((void**)&(c_atom.x), sizeof(MD_FLOAT) * atom->Nmax * 3) ); checkCUDAError( "c_atom.fx malloc", cudaMalloc((void**)&(c_atom.fx), sizeof(MD_FLOAT) * Nlocal) ); @@ -170,25 +179,31 @@ double computeForce( checkCUDAError( "c_atom.sigma6 malloc", cudaMalloc((void**)&(c_atom.sigma6), sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes) ); checkCUDAError( "c_atom.cutforcesq malloc", cudaMalloc((void**)&(c_atom.cutforcesq), sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes) ); - initialized = true; + checkCUDAError( "c_neighs malloc", cudaMalloc((void**)&c_neighs, sizeof(int) * Nlocal * neighbor->maxneighs) ); + checkCUDAError( "c_neigh_numneigh malloc", cudaMalloc((void**)&c_neigh_numneigh, sizeof(int) * Nlocal) ); } - checkCUDAError( "c_atom.x memcpy", cudaMemcpy(c_atom.x, atom->x, sizeof(MD_FLOAT) * atom->Nmax * 3, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.fx memcpy", cudaMemcpy(c_atom.fx, fx, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.fy memcpy", cudaMemcpy(c_atom.fy, fy, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.fz memcpy", cudaMemcpy(c_atom.fz, fz, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.type memcpy", cudaMemcpy(c_atom.type, atom->type, sizeof(int) * atom->Nmax, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.epsilon memcpy", cudaMemcpy(c_atom.epsilon, atom->epsilon, sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.sigma6 memcpy", cudaMemcpy(c_atom.sigma6, atom->sigma6, sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes, cudaMemcpyHostToDevice) ); - checkCUDAError( "c_atom.cutforcesq memcpy", cudaMemcpy(c_atom.cutforcesq, atom->cutforcesq, sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes, cudaMemcpyHostToDevice) ); + printf("2\r\n"); - int *c_neighs; - checkCUDAError( "c_neighs malloc", cudaMalloc((void**)&c_neighs, sizeof(int) * Nlocal * neighbor->maxneighs) ); - checkCUDAError( "c_neighs memcpy", cudaMemcpy(c_neighs, neighbor->neighbors, sizeof(int) * Nlocal * neighbor->maxneighs, cudaMemcpyHostToDevice) ); + if(reneighbourHappenend || !initialized) { + checkCUDAError( "c_atom.x memcpy", cudaMemcpy(c_atom.x, atom->x, sizeof(MD_FLOAT) * atom->Nmax * 3, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.fx memcpy", cudaMemcpy(c_atom.fx, fx, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.fy memcpy", cudaMemcpy(c_atom.fy, fy, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.fz memcpy", cudaMemcpy(c_atom.fz, fz, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.type memcpy", cudaMemcpy(c_atom.type, atom->type, sizeof(int) * atom->Nmax, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.epsilon memcpy", cudaMemcpy(c_atom.epsilon, atom->epsilon, sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.sigma6 memcpy", cudaMemcpy(c_atom.sigma6, atom->sigma6, sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_atom.cutforcesq memcpy", cudaMemcpy(c_atom.cutforcesq, atom->cutforcesq, sizeof(MD_FLOAT) * atom->ntypes * atom->ntypes, cudaMemcpyHostToDevice) ); - int *c_neigh_numneigh; - checkCUDAError( "c_neigh_numneigh malloc", cudaMalloc((void**)&c_neigh_numneigh, sizeof(int) * Nlocal) ); - checkCUDAError( "c_neigh_numneigh memcpy", cudaMemcpy(c_neigh_numneigh, neighbor->numneigh, sizeof(int) * Nlocal, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_neigh_numneigh memcpy", cudaMemcpy(c_neigh_numneigh, neighbor->numneigh, sizeof(int) * Nlocal, cudaMemcpyHostToDevice) ); + checkCUDAError( "c_neighs memcpy", cudaMemcpy(c_neighs, neighbor->neighbors, sizeof(int) * Nlocal * neighbor->maxneighs, cudaMemcpyHostToDevice) ); + } + + printf("3\r\n"); + + printf("4\r\n"); + + printf("5\r\n"); const int num_threads_per_block = num_threads; // this should be multiple of 32 as operations are performed at the level of warps const int num_blocks = ceil((float)Nlocal / (float)num_threads_per_block); @@ -201,10 +216,16 @@ double computeForce( checkCUDAError( "PeekAtLastError", cudaPeekAtLastError() ); checkCUDAError( "DeviceSync", cudaDeviceSynchronize() ); + printf("6\r\n"); + // copy results in c_atom.fx/fy/fz to atom->fx/fy/fz - cudaMemcpy(atom->fx, c_atom.fx, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyDeviceToHost); - cudaMemcpy(atom->fy, c_atom.fy, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyDeviceToHost); - cudaMemcpy(atom->fz, c_atom.fz, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyDeviceToHost); + if(reneighbourHappenend) { + cudaMemcpy(atom->fx, c_atom.fx, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyDeviceToHost); + cudaMemcpy(atom->fy, c_atom.fy, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyDeviceToHost); + cudaMemcpy(atom->fz, c_atom.fz, sizeof(MD_FLOAT) * Nlocal, cudaMemcpyDeviceToHost); + } + + printf("7\r\n"); /* cudaFree(c_atom.x); @@ -215,13 +236,15 @@ double computeForce( cudaFree(c_atom.cutforcesq); */ - cudaFree(c_neighs); cudaFree(c_neigh_numneigh); + // cudaFree(c_neighs); cudaFree(c_neigh_numneigh); cudaProfilerStop(); LIKWID_MARKER_STOP("force"); double E = getTimeStamp(); + initialized = true; + return E-S; } } \ No newline at end of file diff --git a/src/main.c b/src/main.c index 9767bd8..1730dc3 100644 --- a/src/main.c +++ b/src/main.c @@ -111,13 +111,20 @@ double reneighbour( { double S, E; + printf("10.1\r\n"); + S = getTimeStamp(); LIKWID_MARKER_START("reneighbour"); + printf("10.2\r\n"); updateAtomsPbc(atom, param); + printf("10.3\r\n"); setupPbc(atom, param); + printf("10.4\r\n"); updatePbc(atom, param); + printf("10.5\r\n"); //sortAtom(atom); buildNeighbor(atom, neighbor); + printf("10.6\r\n"); LIKWID_MARKER_STOP("reneighbour"); E = getTimeStamp(); @@ -279,6 +286,11 @@ int main(int argc, char** argv) initialIntegrate(¶m, &atom); const bool doReneighbour = (n + 1) % param.every == 0; + const bool doesReneighbourNextRound = (n + 2) % param.every == 0; + + printf("Run %d does reneighbour: %d\r\n", n, doReneighbour); + + printf("10\r\n"); if(doReneighbour) { timer[NEIGH] += reneighbour(¶m, &atom, &neighbor); @@ -286,6 +298,8 @@ int main(int argc, char** argv) updatePbc(&atom, ¶m); } + printf("11\r\n"); + if(param.force_field == FF_EAM) { timer[FORCE] += computeForceEam(&eam, ¶m, &atom, &neighbor, &stats, 0, n + 1); } else { diff --git a/src/neighbor.c b/src/neighbor.c index fbf582d..0abbb26 100644 --- a/src/neighbor.c +++ b/src/neighbor.c @@ -172,6 +172,8 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) { int nall = atom->Nlocal + atom->Nghost; + printf("nall: %d, nmax: %d\r\n", nall, nmax); + /* extend atom arrays if necessary */ if(nall > nmax) { nmax = nall; @@ -183,10 +185,14 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) // neighbor->neighbors = (int*) malloc(nmax * neighbor->maxneighs * sizeof(int*)); } + printf("10.5.1\r\n"); + /* bin local & ghost atoms */ binatoms(atom); int resize = 1; + printf("10.5.2\r\n"); + /* loop over each atom, storing neighbors */ while(resize) { int new_maxneighs = neighbor->maxneighs; @@ -224,6 +230,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) #else const MD_FLOAT cutoff = cutneighsq; #endif + if( rsq <= cutoff ) { neighptr[n++] = j; } @@ -315,9 +322,14 @@ int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin) void binatoms(Atom *atom) { + printf("10.5.1.1\r\n"); int nall = atom->Nlocal + atom->Nghost; int resize = 1; + printf("10.5.1.2\r\n"); + + printf("nall: %d, atom->Nmax: %d\r\n", nall, atom->Nmax); + while(resize > 0) { resize = 0; @@ -325,8 +337,13 @@ void binatoms(Atom *atom) bincount[i] = 0; } + printf("10.5.1.3\r\n"); + for(int i = 0; i < nall; i++) { - int ibin = coord2bin(atom_x(i), atom_y(i), atom_z(i)); + MD_FLOAT x = atom_x(i); + MD_FLOAT y = atom_y(i); + MD_FLOAT z = atom_z(i); + int ibin = coord2bin(x, y, z); if(bincount[ibin] < atoms_per_bin) { int ac = bincount[ibin]++; @@ -336,11 +353,15 @@ void binatoms(Atom *atom) } } + printf("10.5.1.4\r\n"); + if(resize) { free(bins); atoms_per_bin *= 2; bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); } + + printf("10.5.1.5\r\n"); } }