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

First pass at distributed memory for grid

Browse Source
This commit is contained in:
Aidan Thompson
2021-07-02 17:47:45 -06:00
parent 39039d261f
commit e17ace385d
3 changed files with 194 additions and 140 deletions
Download Patch File Download Diff File Expand all files Collapse all files

93
src/SNAP/compute_sna_grid.cpp
View File

@ -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
------------------------------------------------------------------------- */

224
src/compute_grid.cpp
View File

@ -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<int> (comm->xsplit[comm->myloc[0]] * nx);
nxhi = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx) - 1;
if (triclinic) domain->lamda2x(x, x);
nylo = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny);
nyhi = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny) - 1;
nzlo = static_cast<int> (comm->zsplit[comm->myloc[2]] * nz);
nzhi = static_cast<int> (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;
}

17
src/compute_grid.h
View File

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