From e17ace385ddc6d00a485c1c8fd04989526462d6f Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 17:47:45 -0600 Subject: [PATCH] First pass at distributed memory for grid --- src/SNAP/compute_sna_grid.cpp | 93 +++++++------- src/compute_grid.cpp | 224 ++++++++++++++++++++-------------- src/compute_grid.h | 17 ++- 3 files changed, 194 insertions(+), 140 deletions(-) diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index a9f3b6c3fb..196819008c 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -182,66 +182,75 @@ void ComputeSNAGrid::compute_array() snaptr->grow_rij(ntotal); - 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]; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + const int igrid = iz*(nx*ny) + iy*nx + ix; + 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 - // typej = types of neighbors of I within cutoff + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff - int ninside = 0; - for (int j = 0; j < ntotal; j++) { + int ninside = 0; + for (int j = 0; j < ntotal; j++) { - // check that j is in compute group + // check that j is in compute group - if (!(mask[j] & groupbit)) continue; + if (!(mask[j] & groupbit)) continue; - const double delx = xtmp - x[j][0]; - const double dely = ytmp - x[j][1]; - const double delz = ztmp - x[j][2]; - const double rsq = delx*delx + dely*dely + delz*delz; - int jtype = type[j]; - if (rsq < cutsq[jtype][jtype] && rsq>1e-20) { - // printf("ninside = %d\n",ninside); - snaptr->rij[ninside][0] = delx; - snaptr->rij[ninside][1] = dely; - snaptr->rij[ninside][2] = delz; - snaptr->inside[ninside] = j; - snaptr->wj[ninside] = wjelem[jtype]; - snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac; - ninside++; - } - } + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type[j]; + if (rsq < cutsq[jtype][jtype] && rsq>1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac; + ninside++; + } + } - snaptr->compute_ui(ninside); - snaptr->compute_zi(); - snaptr->compute_bi(); - for (int icoeff = 0; icoeff < ncoeff; 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; - for (int icoeff = 0; icoeff < ncoeff; icoeff++) { - double bi = snaptr->blist[icoeff]; + snaptr->compute_ui(ninside); + snaptr->compute_zi(); + snaptr->compute_bi(); + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + gridlocal[size_array_cols_base+icoeff][iz][iy][ix] = snaptr->blist[icoeff]; + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bi = snaptr->blist[icoeff]; - // diagonal element of quadratic matrix + // diagonal element of quadratic matrix - grid[igrid][size_array_cols_base+ncount++] = 0.5*bi*bi; + gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = 0.5*bi*bi; // upper-triangular elements of quadratic matrix for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) - grid[igrid][size_array_cols_base+ncount++] = bi*snaptr->blist[jcoeff]; + gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = bi*snaptr->blist[jcoeff]; } } } + + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + const int igrid = iz*(nx*ny) + iy*nx + ix; + for (int j = 0; j < nvalues; j++) + grid[igrid][size_array_cols_base + j] = gridlocal[size_array_cols_base + j][iz][iy][ix]; + } MPI_Allreduce(&grid[0][0],&gridall[0][0],ngrid*size_array_cols,MPI_DOUBLE,MPI_SUM,world); + } + /* ---------------------------------------------------------------------- memory usage ------------------------------------------------------------------------- */ diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index 011dafa030..b16a51b739 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -21,13 +21,14 @@ #include "force.h" #include "memory.h" #include "error.h" +#include "comm.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), grid(NULL), grid_local(NULL) + Compute(lmp, narg, arg), grid(NULL), local_flags(NULL), gridlocal(NULL) { if (narg < 6) error->all(FLERR,"Illegal compute grid command"); @@ -59,7 +60,8 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : ComputeGrid::~ComputeGrid() { memory->destroy(grid); - memory->destroy(grid_local); + memory->destroy(local_flags); + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); } /* ---------------------------------------------------------------------- */ @@ -72,7 +74,118 @@ void ComputeGrid::init() void ComputeGrid::setup() { - + + set_grid_global(); + set_grid_local(); + allocate(); + assign_coords(); + assign_local_flags(); +} + +/* ---------------------------------------------------------------------- + convert global array index to box coords +------------------------------------------------------------------------- */ + +void ComputeGrid::grid2x(int igrid, double *x) +{ + 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; + + if (triclinic) domain->lamda2x(x, x); + +} + +/* ---------------------------------------------------------------------- + check if grid point is local +------------------------------------------------------------------------- */ + +int ComputeGrid::check_local(int igrid) +{ + 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_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]; + } +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_local_flags() +{ + double x[3]; + for (int igrid = 0; igrid < ngrid; igrid++) { + if (check_local(igrid)) + local_flags[igrid] = 1; + else { + local_flags[igrid] = 0; + memset(grid[igrid],0,size_array_cols); + } + } +} + +/* ---------------------------------------------------------------------- + free and reallocate arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::allocate() +{ + // allocate arrays + + memory->destroy(grid); + memory->destroy(local_flags); + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); + memory->create(grid,size_array_rows,size_array_cols,"grid:grid"); + memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall"); + memory->create(local_flags,size_array_rows,"grid:local_flags"); + memory->create4d_offset(gridlocal,size_array_cols,nzlo,nzhi,nylo,nyhi, + nxlo,nxhi,"grid:gridlocal"); + array = gridall; +} + + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void ComputeGrid::set_grid_global() +{ // calculate grid layout triclinic = domain->triclinic; @@ -100,105 +213,31 @@ void ComputeGrid::setup() delx = 1.0/delxinv; dely = 1.0/delyinv; delz = 1.0/delzinv; - - allocate(); - assign_grid_coords(); - assign_grid_local(); } /* ---------------------------------------------------------------------- - convert global array index to box coords + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) ------------------------------------------------------------------------- */ -void ComputeGrid::grid2x(int igrid, double *x) +void ComputeGrid::set_grid_local() { - int iz = igrid / (nx*ny); - igrid -= iz * (nx*ny); - int iy = igrid / nx; - igrid -= iy * nx; - int ix = igrid; + // global indices of grid range from 0 to N-1 + // nlo,nhi = lower/upper limits of the 3d sub-brick of + // global grid that I own without ghost cells - x[0] = ix*delx; - x[1] = iy*dely; - x[2] = iz*delz; + nxlo = static_cast (comm->xsplit[comm->myloc[0]] * nx); + nxhi = static_cast (comm->xsplit[comm->myloc[0]+1] * nx) - 1; - if (triclinic) domain->lamda2x(x, x); + nylo = static_cast (comm->ysplit[comm->myloc[1]] * ny); + nyhi = static_cast (comm->ysplit[comm->myloc[1]+1] * ny) - 1; + nzlo = static_cast (comm->zsplit[comm->myloc[2]] * nz); + nzhi = static_cast (comm->zsplit[comm->myloc[2]+1] * nz) - 1; + + ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); } -/* ---------------------------------------------------------------------- - check if grid point is local -------------------------------------------------------------------------- */ - -int ComputeGrid::check_grid_local(int igrid) -{ - 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]; - } -} - -/* ---------------------------------------------------------------------- - copy coords to global array -------------------------------------------------------------------------- */ - -void ComputeGrid::assign_grid_local() -{ - double x[3]; - 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); - } - } -} - -/* ---------------------------------------------------------------------- - free and reallocate arrays -------------------------------------------------------------------------- */ - -void ComputeGrid::allocate() -{ - // grow global array if necessary - - 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 ------------------------------------------------------------------------- */ @@ -207,6 +246,7 @@ double ComputeGrid::memory_usage() { double nbytes = size_array_rows*size_array_cols * sizeof(double); // grid - nbytes += size_array_rows*sizeof(int); // grid_local + nbytes += size_array_rows*sizeof(int); // local_flags + nbytes += size_array_cols*ngridlocal*sizeof(double); // gridlocal return nbytes; } diff --git a/src/compute_grid.h b/src/compute_grid.h index a1c938db2e..7e1a4d706b 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -30,11 +30,14 @@ class ComputeGrid : public Compute { double memory_usage(); protected: - int nx, ny, nz; // grid dimensions - int ngrid; // number of grid points + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int ngrid; // number of global grid points + int ngridlocal; // number of local grid points int nvalues; // number of values per grid point double **grid; // global grid double **gridall; // global grid summed over procs + double ****gridlocal; // local grid 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) @@ -43,12 +46,14 @@ class ComputeGrid : public Compute { int nargbase; // number of base class args double cutmax; // largest cutoff distance int size_array_cols_base; // number of columns used for coords, etc. - int *grid_local; // local flag for each grid point + int *local_flags; // 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 + void assign_coords(); // assign coords for grid + void assign_local_flags(); // set local flag for each grid point + int check_local(int); // check if grid point is local + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid private: };