ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Completed brute force parallel implementation using MPI_Allreduce()

Browse Source
This commit is contained in:
Aidan Thompson
2019-10-24 19:48:41 -06:00
parent b44b1f94b7
commit 2fa9e5fefb
4 changed files with 93 additions and 217 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
src/SNAP/compute_sna_grid.cpp
View File

@ -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;
}

4
src/SNAP/compute_sna_grid.h
View File

@ -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;

220
src/compute_grid.cpp
View File

@ -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;
}

26
src/compute_grid.h
View File

@ -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:
};