From 2fa9e5fefb26f3e4540ed7dc6e4e1563e49432bc Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 24 Oct 2019 19:48:41 -0600 Subject: [PATCH] Completed brute force parallel implementation using MPI_Allreduce() --- src/SNAP/compute_sna_grid.cpp | 60 +++------- src/SNAP/compute_sna_grid.h | 4 - src/compute_grid.cpp | 220 ++++++++++------------------------ src/compute_grid.h | 26 ++-- 4 files changed, 93 insertions(+), 217 deletions(-) diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index d9912cdc47..a9f3b6c3fb 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : - ComputeGrid(lmp, narg, arg), cutsq(NULL), list(NULL), sna(NULL), + ComputeGrid(lmp, narg, arg), cutsq(NULL), sna(NULL), radelem(NULL), wjelem(NULL) { double rmin0, rfac0; @@ -120,8 +120,9 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : rmin0,switchflag,bzeroflag); ncoeff = snaptr->ncoeff; - size_array_cols = size_array_cols_base + ncoeff; - if (quadraticflag) size_array_cols += (ncoeff*(ncoeff+1))/2; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2; + size_array_cols = size_array_cols_base + nvalues; array_flag = 1; } @@ -165,26 +166,10 @@ void ComputeSNAGrid::init() /* ---------------------------------------------------------------------- */ -void ComputeSNAGrid::init_list(int /*id*/, NeighList *ptr) -{ - list = ptr; -} - -/* ---------------------------------------------------------------------- */ - void ComputeSNAGrid::compute_array() { invoked_array = update->ntimestep; -// // invoke full neighbor list (will copy or build if necessary) - -// neighbor->build_one(list); - -// const int inum = list->inum; -// const int* const ilist = list->ilist; -// const int* const numneigh = list->numneigh; -// int** const firstneigh = list->firstneigh; - int * const type = atom->type; // compute sna for each gridpoint @@ -197,13 +182,12 @@ void ComputeSNAGrid::compute_array() snaptr->grow_rij(ntotal); - printf("ngridfull = %d\n",ngridfull); - for (int igrid = 0; igrid < ngridfull; igrid++) { - double rtmp[3]; - igridfull2x(igrid, rtmp); - const double xtmp = rtmp[0]; - const double ytmp = rtmp[1]; - const double ztmp = rtmp[2]; + printf("ngrid = %d\n",ngrid); + for (int igrid = 0; igrid < ngrid; igrid++) { + if (!grid_local[igrid]) continue; + const double xtmp = grid[igrid][0]; + const double ytmp = grid[igrid][1]; + const double ztmp = grid[igrid][2]; // rij[][3] = displacements between atom I and those neighbors // inside = indices of neighbors of I within cutoff @@ -237,7 +221,7 @@ void ComputeSNAGrid::compute_array() snaptr->compute_zi(); snaptr->compute_bi(); for (int icoeff = 0; icoeff < ncoeff; icoeff++) - sna[igrid][size_array_cols_base+icoeff] = snaptr->blist[icoeff]; + grid[igrid][size_array_cols_base+icoeff] = snaptr->blist[icoeff]; // printf("igrid = %d %g %g %g %d B0 = %g\n",igrid,xtmp,ytmp,ztmp,ninside,sna[igrid][size_array_cols_base+0]); if (quadraticflag) { int ncount = ncoeff; @@ -246,27 +230,16 @@ void ComputeSNAGrid::compute_array() // diagonal element of quadratic matrix - sna[igrid][ncount++] = 0.5*bi*bi; + grid[igrid][size_array_cols_base+ncount++] = 0.5*bi*bi; // upper-triangular elements of quadratic matrix for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) - sna[igrid][ncount++] = bi*snaptr->blist[jcoeff]; + grid[igrid][size_array_cols_base+ncount++] = bi*snaptr->blist[jcoeff]; } } } - // gather_global_array(); - copy_local_grid(); -} - -/* ---------------------------------------------------------------------- - allocate array in base class and then set up pointers -------------------------------------------------------------------------- */ - -void ComputeSNAGrid::allocate() -{ - ComputeGrid::allocate(); - sna = gridfull; + MPI_Allreduce(&grid[0][0],&gridall[0][0],ngrid*size_array_cols,MPI_DOUBLE,MPI_SUM,world); } /* ---------------------------------------------------------------------- @@ -275,9 +248,8 @@ void ComputeSNAGrid::allocate() double ComputeSNAGrid::memory_usage() { - double bytes = size_array_rows*size_array_cols * sizeof(double); // grid - bytes += snaptr->memory_usage(); // SNA object + double nbytes = snaptr->memory_usage(); // SNA object - return bytes; + return nbytes; } diff --git a/src/SNAP/compute_sna_grid.h b/src/SNAP/compute_sna_grid.h index 2bd18a915e..0242a2962b 100644 --- a/src/SNAP/compute_sna_grid.h +++ b/src/SNAP/compute_sna_grid.h @@ -29,16 +29,12 @@ class ComputeSNAGrid : public ComputeGrid { ComputeSNAGrid(class LAMMPS *, int, char **); ~ComputeSNAGrid(); void init(); - void init_list(int, class NeighList *); void compute_array(); double memory_usage(); - void allocate(); private: - int nmax; int ncoeff; double **cutsq; - class NeighList *list; double **sna; double rcutfac; double *radelem; diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index 63e4a8c885..2f67d22182 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -50,7 +50,7 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : nargbase = iarg - iarg0; - size_array_rows = nx*ny*nz; + size_array_rows = ngrid = nx*ny*nz; size_array_cols_base = 3; } @@ -58,6 +58,8 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : ComputeGrid::~ComputeGrid() { + memory->destroy(grid); + memory->destroy(grid_local); } /* ---------------------------------------------------------------------- */ @@ -78,9 +80,13 @@ void ComputeGrid::setup() if (triclinic == 0) { prd = domain->prd; boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; } else { prd = domain->prd_lamda; boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; } double xprd = prd[0]; @@ -95,78 +101,23 @@ void ComputeGrid::setup() dely = 1.0/delyinv; delz = 1.0/delzinv; - // sufficient conditions for stencil bounding rcut - - // require |delz*mz|^2 <= rcut^2 - // require |dely*my|^2 <= rcut^2 + |delyz*mz_max|^2 - // require |delx*mx|^2 <= rcut^2 + |delxz*mz_max|^2 + |delxy*my_max|^2 - - double delxy = domain->xy/ny; - double delxz = domain->xz/nz; - double delyz = domain->yz/nz; - - if (!triclinic) { - mz = cutmax*delzinv + 1; - my = sqrt(cutmax*cutmax + pow(delyz*mz,2))*delyinv + 1; - mx = sqrt(cutmax*cutmax + pow(delxz*mz,2) - + pow(delxy*my,2))*delxinv + 1; - } else { - double delxinvtmp = nx/domain->xprd; - double delyinvtmp = ny/domain->yprd; - double delzinvtmp = nz/domain->zprd; - mz = cutmax*delzinvtmp + 1; - my = sqrt(cutmax*cutmax + pow(delyz*mz,2))*delyinvtmp + 1; - mx = sqrt(cutmax*cutmax + pow(delxz*mz,2) - + pow(delxy*my,2))*delxinvtmp + 1; - } - - // size global grid to accomodate periodic interactions - - nxfull = nx + 2*mx; - nyfull = ny + 2*my; - nzfull = nz + 2*mz; - nxyfull = nxfull * nyfull; - - x0full = boxlo[0] - mx*delx; - y0full = boxlo[1] - my*dely; - z0full = boxlo[2] - mz*delz; - allocate(); - assign_coords_array(); + assign_grid_coords(); + assign_grid_local(); } /* ---------------------------------------------------------------------- - convert grid index to box coords + convert global array index to box coords ------------------------------------------------------------------------- */ -void ComputeGrid::igridfull2x(int igrid, double *x) +void ComputeGrid::grid2x(int igrid, double *x) { - int iz = igrid / nxyfull; - igrid -= iz * nxyfull; - int iy = igrid / nxfull; - igrid -= iy * nxfull; + int iz = igrid / (nx*ny); + igrid -= iz * (nx*ny); + int iy = igrid / nx; + igrid -= iy * nx; int ix = igrid; - x[0] = x0full+ix*delx; - x[1] = y0full+iy*dely; - x[2] = z0full+iz*delz; - - if (triclinic) domain->lamda2x(x, x); - -} - -/* ---------------------------------------------------------------------- - convert array index to box coords -------------------------------------------------------------------------- */ - -void ComputeGrid::iarray2x(int iarray, double *x) -{ - int iz = iarray / (nx*ny); - iarray -= iz * (nx*ny); - int iy = iarray / nx; - iarray -= iy * nx; - int ix = iarray; - x[0] = ix*delx; x[1] = iy*dely; x[2] = iz*delz; @@ -176,16 +127,43 @@ void ComputeGrid::iarray2x(int iarray, double *x) } /* ---------------------------------------------------------------------- - copy local grid to global array + check if grid point is local ------------------------------------------------------------------------- */ -void ComputeGrid::copy_local_grid() +int ComputeGrid::check_grid_local(int igrid) { - int igridfull; - for (int iarray = 0; iarray < size_array_rows; iarray++) { - igridfull = iarray2igridfull(iarray); - for (int icol = size_array_cols_base; icol < size_array_cols; icol++) - array[iarray][icol] = gridfull[igridfull][icol]; + double x[3]; + + int iz = igrid / (nx*ny); + igrid -= iz * (nx*ny); + int iy = igrid / nx; + igrid -= iy * nx; + int ix = igrid; + + x[0] = ix*delx; + x[1] = iy*dely; + x[2] = iz*delz; + + int islocal = + x[0] >= sublo[0] && x[0] < subhi[0] && + x[1] >= sublo[1] && x[1] < subhi[1] && + x[2] >= sublo[2] && x[2] < subhi[2]; + + return islocal; +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_grid_coords() +{ + double x[3]; + for (int igrid = 0; igrid < ngrid; igrid++) { + grid2x(igrid,x); + grid[igrid][0] = x[0]; + grid[igrid][1] = x[1]; + grid[igrid][2] = x[2]; } } @@ -193,101 +171,33 @@ void ComputeGrid::copy_local_grid() copy coords to global array ------------------------------------------------------------------------- */ -void ComputeGrid::assign_coords_array() +void ComputeGrid::assign_grid_local() { double x[3]; - for (int iarray = 0; iarray < size_array_rows; iarray++) { - iarray2x(iarray,x); - array[iarray][0] = x[0]; - array[iarray][1] = x[1]; - array[iarray][2] = x[2]; + for (int igrid = 0; igrid < ngrid; igrid++) { + if (check_grid_local(igrid)) + grid_local[igrid] = 1; + else { + grid_local[igrid] = 0; + memset(grid[igrid],0,size_array_cols); + } } } -/* ---------------------------------------------------------------------- - gather global array from full grid -------------------------------------------------------------------------- */ - -void ComputeGrid::gather_global_array() -{ - int iarray; - memset(&array[0][0],0,size_array_rows*size_array_cols*sizeof(double)); - - for (int igrid = 0; igrid < ngridfull; igrid++) { - - // inefficient, should exploit shared ix structure - - iarray = igridfull2iarray(igrid); - for (int icol = size_array_cols_base; icol < size_array_cols; icol++) - array[iarray][icol] += gridfull[igrid][icol]; - } -} - -/* ---------------------------------------------------------------------- - convert full grid index to compute array index - inefficient, should exploit shared ix structure -------------------------------------------------------------------------- */ - -int ComputeGrid::igridfull2iarray(int igrid) -{ - int iz = igrid / nxyfull; - igrid -= iz*nxyfull; - int iy = igrid / nxfull; - igrid -= iy*nxfull; - int ix = igrid; - - ix -= mx; - iy -= my; - iz -= mz; - - while (ix < 0) ix += nx; - while (iy < 0) iy += ny; - while (iz < 0) iz += nz; - - while (ix >= nx) ix -= nx; - while (iy >= ny) iy -= ny; - while (iz >= nz) iz -= nz; - - int iarray = (iz * ny + iy) * nx + ix; - - return iarray; -} - -/* ---------------------------------------------------------------------- - convert compute array index to full grid index - inefficient, should exploit shared ix structure -------------------------------------------------------------------------- */ - -int ComputeGrid::iarray2igridfull(int iarray) -{ - int iz = iarray / (nx*ny); - iarray -= iz*(nx*ny); - int iy = iarray / nx; - iarray -= iy*nx; - int ix = iarray; - - ix += mx; - iy += my; - iz += mz; - - int igrid = (iz * nyfull + iy) * nxfull + ix; - - return igrid; -} - /* ---------------------------------------------------------------------- free and reallocate arrays ------------------------------------------------------------------------- */ void ComputeGrid::allocate() { - ngridfull = nxfull*nyfull*nzfull; - // grow global array if necessary - memory->destroy(array); - memory->create(array,size_array_rows,size_array_cols,"sna/grid:array"); - memory->create(gridfull,ngridfull,size_array_cols,"sna/grid:gridfull"); + memory->destroy(grid); + memory->destroy(grid_local); + memory->create(grid,size_array_rows,size_array_cols,"grid:grid"); + memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall"); + memory->create(grid_local,size_array_rows,"grid:grid_local"); + array = gridall; } /* ---------------------------------------------------------------------- memory usage of local data @@ -295,6 +205,8 @@ void ComputeGrid::allocate() double ComputeGrid::memory_usage() { - int nbytes = 0; + double nbytes = size_array_rows*size_array_cols * + sizeof(double); // grid + nbytes += size_array_rows*sizeof(int); // grid_local return nbytes; } diff --git a/src/compute_grid.h b/src/compute_grid.h index 1468f9ad91..a1c938db2e 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -31,28 +31,24 @@ class ComputeGrid : public Compute { protected: int nx, ny, nz; // grid dimensions - int nxfull, nyfull, nzfull; // grid dimensions with ghost points - int nxyfull; // nx_full*ny_full - int ngridfull; // number of full grid points - double **gridfull; // full grid points - int mx, my, mz; // cutmax stencil dimensions + int ngrid; // number of grid points + int nvalues; // number of values per grid point + double **grid; // global grid + double **gridall; // global grid summed over procs int triclinic; // triclinic flag double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri) double delxinv,delyinv,delzinv; // inverse grid spacing double delx,dely,delz; // grid spacing - double x0full, y0full, z0full; // origin of full grid int nargbase; // number of base class args double cutmax; // largest cutoff distance int size_array_cols_base; // number of columns used for coords, etc. - virtual void allocate(); - void igridfull2x(int, double*); // convert full grid point to coord - void iarray2x(int, double*); // convert array point to coord - void gather_global_array(); // gather global array from full grid - void copy_local_grid(); // copy local grid to global array - void assign_coords_array(); // assign coords to global array - int igridfull2iarray(int); // convert full grid index to compute array index - int iarray2igridfull(int); // convert compute array index to full grid index - + int *grid_local; // local flag for each grid point + void allocate(); + void grid2x(int, double*); // convert grid point to coord + void assign_grid_coords(); // assign coords for grid + void assign_grid_local(); // set local flag for each grid point + int check_grid_local(int); // check if grid point is local private: };