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

support for more caller options in Grid2d/3d

Browse Source
This commit is contained in:
Steve Plimpton
2022-11-02 15:01:58 -06:00
parent 7346aee4ad
commit 8af384243f
5 changed files with 357 additions and 325 deletions
Download Patch File Download Diff File Expand all files Collapse all files

81
src/KSPACE/msm.cpp
View File

@ -594,32 +594,17 @@ void MSM::allocate()
memory->create2d_offset(phi1d,3,-order,order,"msm:phi1d");
memory->create2d_offset(dphi1d,3,-order,order,"msm:dphi1d");
// one Grid3d using all processors for finest grid level
gcall = new Grid3d(lmp,world,nx_msm[0],ny_msm[0],nz_msm[0]);
gcall->set_distance(0.5*neighbor->skin);
gcall->set_stencil_atom(-nlower,nupper);
gcall->setup_grid(nxlo_in[0],nxhi_in[0],nylo_in[0],
nyhi_in[0],nzlo_in[0],nzhi_in[0],
nxlo_out_all,nxhi_out_all,nylo_out_all,
nyhi_out_all,nzlo_out_all,nzhi_out_all);
gcall->set_larger_grid(nxlo_out[0],nxhi_out[0],nylo_out[0],
nyhi_out[0],nzlo_out[0],nzhi_out[0]);
// one Grid3d for finest grid level, using world comm and all procs
// use set_caller_grid() b/c MSM allocates local grid > out_all values
// NOTE: or is it out[0] ??
// NOTE: worry about flag = 1 here, 2 later
/*
gcall = new Grid3d(lmp,world,1,nx_msm[0],ny_msm[0],nz_msm[0],
gcall = new Grid3d(lmp,world,nx_msm[0],ny_msm[0],nz_msm[0],
nxlo_in[0],nxhi_in[0],nylo_in[0],
nyhi_in[0],nzlo_in[0],nzhi_in[0],
nxlo_out_all,nxhi_out_all,nylo_out_all,
nyhi_out_all,nzlo_out_all,nzhi_out_all,
nxlo_out[0],nxhi_out[0],nylo_out[0],
nyhi_out[0],nzlo_out[0],nzhi_out[0]);
*/
nyhi_out_all,nzlo_out_all,nzhi_out_all);
gcall->set_caller_grid(nxlo_out[0],nxhi_out[0],nylo_out[0],
nyhi_out[0],nzlo_out[0],nzhi_out[0]);
gcall->setup_comm(ngcall_buf1,ngcall_buf2);
npergrid = 1;
@ -639,34 +624,20 @@ void MSM::allocate()
memory->create3d_offset(egrid[n],nzlo_out[n],nzhi_out[n],
nylo_out[n],nyhi_out[n],nxlo_out[n],nxhi_out[n],"msm:egrid");
// create commgrid object for rho and electric field communication
// one Grid3d per level, using level-specific comm for coarser grids
// use set_proc_neigh() b/c MSM excludes non-participating procs from owned/ghost comm
if (active_flag[n]) {
delete gc[n];
// NOTE: why is n = 0 same as all for grid size ?
gc[n] = new Grid3d(lmp,world_levels[n],nx_msm[n],ny_msm[n],nz_msm[n]);
gc[n]->set_stencil_atom(-nlower,nupper);
gc[n]->setup_grid(nxlo_in[n],nxhi_in[n],nylo_in[n],
nyhi_in[n],nzlo_in[n],nzhi_in[n],
nxlo_out[n],nxhi_out[n],nylo_out[n],
nyhi_out[n],nzlo_out[n],nzhi_out[n]);
gc[n] = new Grid3d(lmp,world_levels[n],nx_msm[n],ny_msm[n],nz_msm[n],
nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n],nzlo_in[n],nzhi_in[n],
nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n],nzlo_out[n],nzhi_out[n]);
int **procneigh = procneigh_levels[n];
gc[n]->set_proc_neighs(procneigh[0][0],procneigh[0][1],procneigh[1][0],
procneigh[1][1],procneigh[2][0],procneigh[2][1]);
/*
gc[n] = new Grid3d(lmp,world_levels[n],2,nx_msm[n],ny_msm[n],nz_msm[n],
nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n],
nzlo_in[n],nzhi_in[n],
nxlo_out[n],nxhi_out[n],nylo_out[n],nyhi_out[n],
procneigh[0][0],procneigh[0][1],procneigh[1][0],
procneigh[1][1],procneigh[2][0],procneigh[2][1]);
*/
gc[n]->setup_comm(ngc_buf1[n],ngc_buf2[n]);
npergrid = 1;
memory->destroy(gc_buf1[n]);
@ -1194,18 +1165,23 @@ void MSM::set_grid_local()
for (int n = 0; n < levels; n++) {
// deleted and nullify grid arrays since the number or offset of gridpoints may change
// delete and nullify grid arrays since the number or offset of gridpoints may change
memory->destroy3d_offset(qgrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
memory->destroy3d_offset(egrid[n],nzlo_out[n],nylo_out[n],nxlo_out[n]);
// partition global grid across procs
// n xyz lo/hi in[] = lower/upper bounds of global grid this proc owns
// indices range from 0 to N-1 inclusive in each dim
// use Grid3d to partition each grid into owned cells on each proc
// assignment of ghost cells is done below
// Grid3d is called later in allocate() with owned+ghost bounds
//comm->partition_grid(nx_msm[n],ny_msm[n],nz_msm[n],0.0,
// nxlo_in[n],nxhi_in[n],nylo_in[n],nyhi_in[n],
// nzlo_in[n],nzhi_in[n]);
gcall = new Grid3d(lmp,world,nx_msm[n],ny_msm[n],nz_msm[n]);
gcall->setup_grid(nxlo_in[n],nxhi_in[n],nylo_in[n],
nyhi_in[n],nzlo_in[n],nzhi_in[n],
nxlo_out[n],nxhi_out[n],nylo_out[n],
nyhi_out[n],nzlo_out[n],nzhi_out[n]);
delete gcall;
// nlower/nupper = stencil size for mapping particles to grid
nlower = -(order-1)/2;
nupper = order/2;
@ -1236,10 +1212,12 @@ void MSM::set_grid_local()
// position a particle in my box can be at
// dist[3] = particle position bound = subbox + skin/2.0
// nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping
// for n = 0, use a smaller ghost region for interpolation
// for n > 0, larger ghost region needed for direct sum and restriction/prolongation
double dist[3];
double cuthalf = 0.0;
if (n == 0) cuthalf = 0.5*neighbor->skin; // only applies to finest grid
if (n == 0) cuthalf = 0.5*neighbor->skin; // only applies to finest grid
dist[0] = dist[1] = dist[2] = cuthalf;
if (triclinic) MathExtra::tribbox(domain->h,cuthalf,&dist[0]);
@ -1249,14 +1227,12 @@ void MSM::set_grid_local()
nx_msm[n]/xprd + OFFSET) - OFFSET;
nhi = static_cast<int> ((subhi[0]+dist[0]-boxlo[0]) *
nx_msm[n]/xprd + OFFSET) - OFFSET;
if (n == 0) {
// use a smaller ghost region for interpolation
nxlo_out_all = nlo + nlower;
nxhi_out_all = nhi + nupper;
}
// a larger ghost region is needed for the direct sum and for restriction/prolongation
nxlo_out[n] = nlo + MIN(-order,nxlo_direct);
nxhi_out[n] = nhi + MAX(order,nxhi_direct);
@ -1268,6 +1244,7 @@ void MSM::set_grid_local()
nylo_out_all = nlo + nlower;
nyhi_out_all = nhi + nupper;
}
nylo_out[n] = nlo + MIN(-order,nylo_direct);
nyhi_out[n] = nhi + MAX(order,nyhi_direct);

344
src/grid2d.cpp
View File

@ -40,9 +40,11 @@ static constexpr int OFFSET = 16384;
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
constructor to create a 3d distributed grid
constructor to create a 2d distributed grid
Grid3d assigns owned/ghost cells to each proc via setup_grid()
it MUST be called after constructor
gcomm = caller's communicator
gnx,gny = global grid size
gnx,gny,gnz = global grid size
------------------------------------------------------------------------- */
Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny) : Pointers(lmp)
@ -55,32 +57,35 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny) : Pointers(lmp)
ny = gny;
// default settings, can be overridden by set() methods
// these affect assignment of owned and ghost cells
maxdist = 0.0;
stencil_grid_lo = stencil_grid_hi = 0;
stencil_atom_lo = stencil_atom_hi = 0;
shift_grid = shift_atom = 0.0;
shift_grid = 0.5;
shift_atom_lo = shift_atom_hi = 0.0;
yextra = 0;
yfactor = 1.0;
}
/* ----------------------------------------------------------------------
constructor called by PPPM classes
gcomm = world communicator
gnx, gny = size of global grid
i xy lohi = portion of global grid this proc owns, 0 <= index < N
o xy lohi = owned grid portion + ghost grid cells needed in all directions
if o indices are < 0 or hi indices are >= N,
then grid is treated as periodic in that dimension,
communication is done across the periodic boundaries
alternate constructor to create a 2d distributed grid
caller assigns owned/ghost cells to each proc
setup_grid() must NOT be called
gcomm = caller's communicator
gnx,gny = global grid size
i xy lo/hi = extent of owned grid cells on this proc
o xy lo/hi = extent of owned+ghost grid cells on this proc
owned and ghost indices are inclusive
owned indices range from 0 to N-1
ghost indices can extend < 0 or >= N
------------------------------------------------------------------------- */
Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm,
int gnx, int gny,
Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny,
int ixlo, int ixhi, int iylo, int iyhi,
int oxlo, int oxhi, int oylo, int oyhi)
: Pointers(lmp)
int oxlo, int oxhi, int oylo, int oyhi) :
Pointers(lmp)
{
// store commnicator and global grid size
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
@ -88,87 +93,22 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm,
nx = gnx;
ny = gny;
// store info about Comm decomposition needed for remap operation
// two Grid instances will exist for duration of remap
// each must know Comm decomp at time Grid instance was created
// store owned/ghost indices provided by caller
inxlo = ixlo;
inxhi = ixhi;
inylo = iylo;
inyhi = iyhi;
extract_comm_info();
outxlo = oxlo;
outxhi = oxhi;
outylo = oylo;
outyhi = oyhi;
// store grid bounds and proc neighs
// additional intialization
// other constructor invokes this from setup_grid()
if (layout != Comm::LAYOUT_TILED) {
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
procxlo,procxhi,procylo,procyhi);
} else {
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
0,0,0,0);
}
}
/* ----------------------------------------------------------------------
constructor called by MSM
gcomm = world communicator or sub-communicator for a hierarchical grid
flag = 1 if e xy lohi values = larger grid stored by caller in gcomm = world
flag = 2 if e xy lohi values = 6 neighbor procs in gcomm
gnx, gny = size of global grid
i xy lohi = portion of global grid this proc owns, 0 <= index < N
o xy lohi = owned grid portion + ghost grid cells needed in all directions
e xy lohi for flag = 1: extent of larger grid stored by caller
e xy lohi for flag = 2: 4 neighbor procs
------------------------------------------------------------------------- */
Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag,
int gnx, int gny,
int ixlo, int ixhi, int iylo, int iyhi,
int oxlo, int oxhi, int oylo, int oyhi,
int exlo, int exhi, int eylo, int eyhi)
: Pointers(lmp)
{
// store commnicator and global grid size
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
// store info about Comm decomposition needed for remap operation
// two Grid instances will exist for duration of remap
// each must know Comm decomp at time Grid instance was created
extract_comm_info();
// store grid bounds and proc neighs
if (flag == 1) {
if (layout != Comm::LAYOUT_TILED) {
// this assumes gcomm = world
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi,
procxlo,procxhi,procylo,procyhi);
} else {
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi,
0,0,0,0);
}
} else if (flag == 2) {
if (layout != Comm::LAYOUT_TILED) {
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi);
} else {
error->all(FLERR,"Grid2d does not support tiled layout for MSM");
}
}
initialize();
}
/* ---------------------------------------------------------------------- */
@ -213,42 +153,13 @@ Grid2d::~Grid2d()
deallocate_remap();
}
/* ----------------------------------------------------------------------
store grid bounds and proc neighs from caller in internal variables
------------------------------------------------------------------------- */
void Grid2d::store(int ixlo, int ixhi, int iylo, int iyhi,
int oxlo, int oxhi, int oylo, int oyhi,
int fxlo, int fxhi, int fylo, int fyhi,
int pxlo, int pxhi, int pylo, int pyhi)
{
inxlo = ixlo;
inxhi = ixhi;
inylo = iylo;
inyhi = iyhi;
outxlo = oxlo;
outxhi = oxhi;
outylo = oylo;
outyhi = oyhi;
fullxlo = fxlo;
fullxhi = fxhi;
fullylo = fylo;
fullyhi = fyhi;
procxlo = pxlo;
procxhi = pxhi;
procylo = pylo;
procyhi = pyhi;
}
// ----------------------------------------------------------------------
// set Grid parameters
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
maxdist = max distance outside proc subdomain a particle can be
used to determine extent of ghost cells
------------------------------------------------------------------------- */
void Grid2d::set_distance(double distance)
@ -258,7 +169,6 @@ void Grid2d::set_distance(double distance)
/* ----------------------------------------------------------------------
# of grid cells beyond an owned grid cell that caller accesses
e.g. for a finite different stencial
can be different in lo vs hi direction
------------------------------------------------------------------------- */
@ -269,8 +179,7 @@ void Grid2d::set_stencil_grid(int lo, int hi)
}
/* ----------------------------------------------------------------------
# of grid cells beyond a particle's grid cell that caller updates
e.g. for smearing a point charge to the grid
# of grid cells beyond a particle's grid cell that caller accesses
can be different in lo vs hi direction
------------------------------------------------------------------------- */
@ -281,9 +190,9 @@ void Grid2d::set_stencil_atom(int lo, int hi)
}
/* ----------------------------------------------------------------------
shift_grid = offset of position of grid point within grid cell
0.5 = cell center, 0.0 = lower-left corner of cell
used to determine which proc owns the grid cell within its subdomain
shift_grid = offset within grid cell of position of grid point
0.5 = cell center (default), 0.0 = lower-left corner of cell
used to decide which proc owns a grid cell (grid pt within subdomain)
------------------------------------------------------------------------- */
void Grid2d::set_shift_grid(double shift)
@ -292,15 +201,58 @@ void Grid2d::set_shift_grid(double shift)
}
/* ----------------------------------------------------------------------
shift_atom = offset of atoms when caller maps them to grid cells
shift_atom = offset added to atoms when caller maps them to grid cells
0.5 = half a grid cell, 0.0 = no offset
used when computing possible ghost extent
PPPM uses 0.5 when order is odd, 0.0 when order is even
used to compute maximum possible ghost extents
use of lo/hi allows max ghost extent on each side to be different
------------------------------------------------------------------------- */
void Grid2d::set_shift_atom(double shift)
void Grid2d::set_shift_atom(double shift_lo, double shift_hi)
{
shift_atom = shift;
shift_atom_lo = shift_lo;
shift_atom_hi = shift_hi;
}
/* ----------------------------------------------------------------------
enable extra grid cells in Y
factor = muliplication factor on box size Y and thus grid size
factor > 1.0 when grid extends beyond Y box size (3.0 = tripled in size)
default zextra = 0, factor = 1.0 (no extra grid cells in Y)
------------------------------------------------------------------------- */
void Grid2d::set_yfactor(double factor)
{
if (factor == 1.0) yextra = 0;
else yextra = 1;
yfactor = factor;
}
/* ----------------------------------------------------------------------
set IDs of proc neighbors used in uniform local owned/ghost comm
called AFTER setup_grid() but BEFORE setup_comm() to override
the processor neighbors stored by extract_comm()
------------------------------------------------------------------------- */
void Grid2d::set_proc_neighs(int pxlo, int pxhi, int pylo, int pyhi)
{
procxlo = pxlo;
procxhi = pxhi;
procylo = pylo;
procyhi = pyhi;
}
/* ----------------------------------------------------------------------
set allocation dimensions of caller grid used by indices() to setup pack/unpack
called AFTER setup_grid() but BEFORE setup_comm() to override
the caller grid size set by setup_grid() and used in indices()
------------------------------------------------------------------------- */
void Grid2d::set_caller_grid(int fxlo, int fxhi, int fylo, int fyhi)
{
fullxlo = fxlo;
fullxhi = fxhi;
fullylo = fylo;
fullyhi = fyhi;
}
// ----------------------------------------------------------------------
@ -382,32 +334,56 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
if (comm->layout != Comm::LAYOUT_TILED) {
fraclo = comm->xsplit[comm->myloc[0]];
frachi = comm->xsplit[comm->myloc[0]+1];
partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi);
partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi);
fraclo = comm->ysplit[comm->myloc[1]];
frachi = comm->ysplit[comm->myloc[1]+1];
partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi);
partition_grid(ny,fraclo,frachi,shift_grid,yextra,inylo,inyhi);
} else {
fraclo = comm->mysplit[0][0];
frachi = comm->mysplit[0][1];
partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi);
partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi);
fraclo = comm->mysplit[1][0];
frachi = comm->mysplit[1][1];
partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi);
partition_grid(ny,fraclo,frachi,shift_grid,yextra,inylo,inyhi);
}
// extend owned grid bounds with ghost grid cells in each direction
ghost_grid();
// additional intialization
// other constructor invokes this directly
initialize();
// return values
ixlo = inxlo;
ixhi = inxhi;
iylo = inylo;
iyhi = inyhi;
oxlo = outxlo;
oxhi = outxhi;
oylo = outylo;
oyhi = outyhi;
}
/* ----------------------------------------------------------------------
additional one-time setup common to both constructors
---------------------------------------------------------------------- */
void Grid2d::initialize()
{
// error check on size of grid stored by this proc
bigint total = (bigint) (outxhi - outxlo + 1) * (outyhi - outylo + 1);
if (total > MAXSMALLINT)
error->one(FLERR, "Too many owned+ghost grid3d points");
error->one(FLERR, "Too many owned+ghost grid2d points");
// default = caller grid is allocated to ghost grid
// used when computing pack/unpack lists in indices()
// NOTE: allow caller to override this
// these values can be overridden using set_caller_grid()
fullxlo = outxlo;
fullxhi = outxhi;
@ -439,18 +415,6 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
// each must know Comm decomp at time Grid instance was created
extract_comm_info();
// return values
ixlo = inxlo;
ixhi = inxhi;
iylo = inylo;
iyhi = inyhi;
oxlo = outxlo;
oxhi = outxhi;
oylo = outylo;
oyhi = outyhi;
}
/* ----------------------------------------------------------------------
@ -460,11 +424,10 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
ngrid = extent of global grid in a dimension
indices into the global grid range from 0 to Ngrid-1 in that dim
shift factor determines position of grid pt within grid cell
// NOTE: for future support of zfactor
zfactor = 0.0 if the grid exactly covers the simulation box
zfactor > 1.0 if the grid extends beyond the +z boundary by this factor
used by 2d slab-mode PPPM
this effectively maps proc sub-grids to a smaller subset of the grid
shift = 0.5 for cell center, 0.0 for lower-left corner
extra = 0 if grid exactly covers the simulation box
extra = 1 if grid extends beyond the +y boundary by yfactor
effectively maps proc partitions to the box-size subset of the grid
lo/hi = inclusive lo/hi bounds for brick of global grid cells I own
lo grid index = first grid pt >= fraclo bound
hi grid index = last grid pt < frachi bound
@ -474,12 +437,19 @@ void Grid2d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
------------------------------------------------------------------------- */
void Grid2d::partition_grid(int ngrid, double fraclo, double frachi,
double shift, int &lo, int &hi)
double shift, int extra, int &lo, int &hi)
{
lo = static_cast<int> (fraclo * ngrid);
while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
hi = static_cast<int> (frachi * ngrid);
while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
if (extra == 0) {
lo = static_cast<int> (fraclo * ngrid);
while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
hi = static_cast<int> (frachi * ngrid);
while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
} else {
lo = static_cast<int> (fraclo * ngrid/yfactor);
while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
hi = static_cast<int> (frachi * ngrid/yfactor);
while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
}
}
/* ----------------------------------------------------------------------
@ -521,19 +491,47 @@ void Grid2d::ghost_grid()
double dxinv = nx / prd[0];
double dyinv = ny / prd[1];
double dyinv_slab = ny / (prd[1] * yfactor);
int lo, hi;
lo = static_cast<int>((sublo[0]-dist[0]-boxlo[0]) * dxinv + shift_atom + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[0]+dist[0]-boxlo[0]) * dxinv + shift_atom + OFFSET) - OFFSET;
lo = static_cast<int>((sublo[0]-dist[0]-boxlo[0]) * dxinv + shift_atom_lo + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[0]+dist[0]-boxlo[0]) * dxinv + shift_atom_hi + OFFSET) - OFFSET;
outxlo = MIN(lo - stencil_atom_lo, inxlo - stencil_grid_lo);
outxhi = MAX(hi + stencil_atom_hi, inxhi + stencil_grid_hi);
lo = static_cast<int>((sublo[1]-dist[1]-boxlo[1]) * dyinv + shift_atom + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[1]+dist[1]-boxlo[1]) * dyinv + shift_atom + OFFSET) - OFFSET;
outylo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo);
outyhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi);
if (yextra == 0) {
lo = static_cast<int>((sublo[1]-dist[1]-boxlo[1]) * dyinv + shift_atom_lo + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[1]+dist[1]-boxlo[1]) * dyinv + shift_atom_hi + OFFSET) - OFFSET;
outylo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo);
outyhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi);
} else {
lo = static_cast<int>((sublo[1]-dist[1]-boxlo[1]) * dyinv_slab + shift_atom_lo + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[1]+dist[1]-boxlo[1]) * dyinv_slab + shift_atom_hi + OFFSET) - OFFSET;
outylo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo);
outyhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi);
}
// if yextra = 1:
// adjust grid boundaries for processors at +y end,
// to include added empty grid cells between periodically repeating slabs
// in this case:
// want grid data forward communicated from +y proc to -y proc, but not vice versa
// want grid data reverse communicated from -y proc to +y proc, but not vice versa
// this is accomplished by inyhi = outyhi on +y end (no ghost cells)
// also insure no other procs use ghost cells beyond +y limit
if (yextra) {
if (layout != Comm::LAYOUT_TILED) {
if (comm->myloc[1] == comm->procgrid[1]-1) inyhi = outyhi = ny - 1;
} else {
if (comm->mysplit[1][1] == 1.0) inyhi = outyhi = ny - 1;
}
outyhi = MIN(outyhi,ny-1);
}
// limit out xyz lo/hi indices to global grid for non-periodic dims
// if yextra = 1, grid is still fully periodic
int *periodicity = domain->periodicity;
@ -558,7 +556,7 @@ void Grid2d::extract_comm_info()
layout = comm->layout;
// for non TILED layout:
// proc xyz lohi = my 6 neighbor procs in this MPI_Comm
// proc xyz lohi = my 64neighbor procs in this MPI_Comm
// NOTE: will need special logic for MSM case with different MPI_Comm
// xyz split = copy of 1d vectors in Comm
// grid2proc = copy of 3d array in Comm
@ -1651,8 +1649,8 @@ int Grid2d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
for (int j = jproclo; j <= jprochi; j++)
for (int i = iproclo; i <= iprochi; i++) {
partition_grid(nx,xsplit[i],xsplit[i+1],shift_grid,obox[0],obox[1]);
partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,obox[2],obox[3]);
partition_grid(nx,xsplit[i],xsplit[i+1],shift_grid,0,obox[0],obox[1]);
partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,yextra,obox[2],obox[3]);
if (noverlap_list == maxoverlap_list) grow_overlap();
overlap_list[noverlap_list].proc = grid2proc[i][j][0];

21
src/grid2d.h
View File

@ -23,18 +23,17 @@ class Grid2d : protected Pointers {
enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes
Grid2d(class LAMMPS *, MPI_Comm, int, int);
Grid2d(class LAMMPS *, MPI_Comm, int, int, int, int, int, int, int, int, int, int);
Grid2d(class LAMMPS *, MPI_Comm, int, int, int, int, int, int, int, int, int, int,
int, int, int, int, int);
~Grid2d() override;
void set_distance(double);
void set_stencil_grid(int, int);
void set_stencil_atom(int, int);
void set_shift_grid(double);
void set_shift_atom(double);
void set_shift_atom(double, double);
void set_yfactor(double);
void set_caller_grid(int, int, int, int);
void set_proc_neighs(int, int, int, int);
int identical(Grid2d *);
void get_size(int &, int &);
@ -68,8 +67,12 @@ protected:
int stencil_grid_lo,stencil_grid_hi; // grid cells accessed beyond owned cell
double shift_grid; // location of grid point within grid cell
// only affects which proc owns grid cell
double shift_atom; // shift applied to atomd when mapped to grid cell by caller
double shift_atom_lo,shift_atom_hi;; // max shift applied to atoms
// when mapped to grid cell by caller
// can be different in lo/hi directions
// only affects extent of ghost cells
int yextra; // 1 if extra grid cells in Y, 0 if not
double yfactor; // multiplier on extent of grid in Y direction
// extent of my owned and ghost cells
@ -224,13 +227,11 @@ protected:
// internal methods
// -------------------------------------------
void partition_grid(int, double, double, double, int &, int &);
void initialize();
void partition_grid(int, double, double, double, int, int &, int &);
void ghost_grid();
void extract_comm_info();
void store(int, int, int, int, int, int, int, int,
int, int, int, int, int, int, int, int);
void setup_comm_brick(int &, int &);
void setup_comm_tiled(int &, int &);
int ghost_adjacent_brick();

224
src/grid3d.cpp
View File

@ -41,6 +41,8 @@ static constexpr int OFFSET = 16384;
/* ----------------------------------------------------------------------
constructor to create a 3d distributed grid
Grid3d assigns owned/ghost cells to each proc via setup_grid()
it MUST be called after constructor
gcomm = caller's communicator
gnx,gny,gnz = global grid size
------------------------------------------------------------------------- */
@ -57,6 +59,7 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) :
nz = gnz;
// default settings, can be overridden by set() methods
// these affect assignment of owned and ghost cells
maxdist = 0.0;
stencil_grid_lo = stencil_grid_hi = 0;
@ -64,6 +67,56 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) :
shift_grid = 0.5;
shift_atom_lo = shift_atom_hi = 0.0;
zextra = 0;
zfactor = 1.0;
}
/* ----------------------------------------------------------------------
alternate constructor to create a 3d distributed grid
caller assigns owned/ghost cells to each proc
setup_grid() must NOT be called
used by MSM b/c its definition of ghost cells is complex
gcomm = caller's communicator
gnx,gny,gnz = global grid size
i xyz lo/hi = extent of owned grid cells on this proc
o xyz lo/hi = extent of owned+ghost grid cells on this proc
owned and ghost indices are inclusive
owned indices range from 0 to N-1
ghost indices can extend < 0 or >= N
------------------------------------------------------------------------- */
Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz,
int ixlo, int ixhi, int iylo, int iyhi, int izlo, int izhi,
int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi) :
Pointers(lmp)
{
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
nz = gnz;
// store owned/ghost indices provided by caller
inxlo = ixlo;
inxhi = ixhi;
inylo = iylo;
inyhi = iyhi;
inzlo = izlo;
inzhi = izhi;
outxlo = oxlo;
outxhi = oxhi;
outylo = oylo;
outyhi = oyhi;
outzlo = ozlo;
outzhi = ozhi;
// additional intialization
// other constructor invokes this from setup_grid()
initialize();
}
/* ---------------------------------------------------------------------- */
@ -150,7 +203,7 @@ void Grid3d::set_stencil_atom(int lo, int hi)
/* ----------------------------------------------------------------------
shift_grid = offset within grid cell of position of grid point
0.5 = cell center (default), 0.0 = lower-left corner of cell
used by Grid to decide which proc owns a grid cell (grid pt within subdomain)
used to decide which proc owns a grid cell (grid pt within subdomain)
------------------------------------------------------------------------- */
void Grid3d::set_shift_grid(double shift)
@ -161,7 +214,7 @@ void Grid3d::set_shift_grid(double shift)
/* ----------------------------------------------------------------------
shift_atom = offset added to atoms when caller maps them to grid cells
0.5 = half a grid cell, 0.0 = no offset
used by Grid to compute maximum possible ghost extents
used to compute maximum possible ghost extents
use of lo/hi allows max ghost extent on each side to be different
PPPM uses 0.5 when stencil order is odd, 0.0 when order is even
PPPM/stagger applies different shift values for 2 stagger iterations
@ -190,7 +243,7 @@ void Grid3d::set_zfactor(double factor)
/* ----------------------------------------------------------------------
set IDs of proc neighbors used in uniform local owned/ghost comm
call this AFTER setup_grid() but BEFORE setup_comm() to override
called AFTER setup_grid() but BEFORE setup_comm() to override
the processor neighbors stored by extract_comm()
used by MSM to exclude non-participating procs for coarse grid comm
------------------------------------------------------------------------- */
@ -208,13 +261,13 @@ void Grid3d::set_proc_neighs(int pxlo, int pxhi, int pylo, int pyhi,
/* ----------------------------------------------------------------------
set allocation dimensions of caller grid used by indices() to setup pack/unpack
call this AFTER setup_grid() but BEFORE setup_comm() to override
the caller grid size set by setup_grid()
called AFTER setup_grid() but BEFORE setup_comm() to override
the caller grid size set by setup_grid() and used in indices()
used by MSM to allow a larger level 0 grid to be allocated
with more ghost cells for other operations
------------------------------------------------------------------------- */
void Grid3d::set_larger_grid(int fxlo, int fxhi, int fylo, int fyhi,
void Grid3d::set_caller_grid(int fxlo, int fxhi, int fylo, int fyhi,
int fzlo, int fzhi)
{
fullxlo = fxlo;
@ -308,36 +361,63 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
int &ozlo, int &ozhi)
{
// owned grid cells = those whose grid point is within proc subdomain
// shift_grid = 0.5 for grid point at cell center, 0.0 for lower-left corner
double fraclo,frachi;
if (comm->layout != Comm::LAYOUT_TILED) {
fraclo = comm->xsplit[comm->myloc[0]];
frachi = comm->xsplit[comm->myloc[0]+1];
partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi);
partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi);
fraclo = comm->ysplit[comm->myloc[1]];
frachi = comm->ysplit[comm->myloc[1]+1];
partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi);
partition_grid(ny,fraclo,frachi,shift_grid,0,inylo,inyhi);
fraclo = comm->zsplit[comm->myloc[2]];
frachi = comm->zsplit[comm->myloc[2]+1];
partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi);
partition_grid(nz,fraclo,frachi,shift_grid,zextra,inzlo,inzhi);
} else {
fraclo = comm->mysplit[0][0];
frachi = comm->mysplit[0][1];
partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi);
partition_grid(nx,fraclo,frachi,shift_grid,0,inxlo,inxhi);
fraclo = comm->mysplit[1][0];
frachi = comm->mysplit[1][1];
partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi);
partition_grid(ny,fraclo,frachi,shift_grid,0,inylo,inyhi);
fraclo = comm->mysplit[2][0];
frachi = comm->mysplit[2][1];
partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi);
partition_grid(nz,fraclo,frachi,shift_grid,zextra,inzlo,inzhi);
}
// extend owned grid bounds with ghost grid cells in each direction
ghost_grid();
// additional intialization
// other constructor invokes this directly
initialize();
// return values
ixlo = inxlo;
ixhi = inxhi;
iylo = inylo;
iyhi = inyhi;
izlo = inzlo;
izhi = inzhi;
oxlo = outxlo;
oxhi = outxhi;
oylo = outylo;
oyhi = outyhi;
ozlo = outzlo;
ozhi = outzhi;
}
/* ----------------------------------------------------------------------
additional one-time setup common to both constructors
---------------------------------------------------------------------- */
void Grid3d::initialize()
{
// error check on size of grid stored by this proc
bigint total = (bigint)
@ -347,7 +427,7 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
// default = caller grid is allocated to ghost grid
// used when computing pack/unpack lists in indices()
// these values can be overridden by caller using set_larger_grid()
// these values can be overridden using set_caller_grid()
fullxlo = outxlo;
fullxhi = outxhi;
@ -381,22 +461,6 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
// each must know Comm decomp at time Grid instance was created
extract_comm_info();
// return values
ixlo = inxlo;
ixhi = inxhi;
iylo = inylo;
iyhi = inyhi;
izlo = inzlo;
izhi = inzhi;
oxlo = outxlo;
oxhi = outxhi;
oylo = outylo;
oyhi = outyhi;
ozlo = outzlo;
ozhi = outzhi;
}
/* ----------------------------------------------------------------------
@ -406,11 +470,10 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
ngrid = extent of global grid in a dimension
indices into the global grid range from 0 to Ngrid-1 in that dim
shift factor determines position of grid pt within grid cell
// NOTE: for future support of zfactor
zfactor = 0.0 if the grid exactly covers the simulation box
zfactor > 1.0 if the grid extends beyond the +z boundary by this factor
used by 2d slab-mode PPPM
this effectively maps proc sub-grids to a smaller subset of the grid
shift = 0.5 for cell center, 0.0 for lower-left corner
extra = 0 if grid exactly covers the simulation box
extra = 1 if grid extends beyond the +z boundary by zfactor (PPPM slab)
effectively maps proc partitions to the box-size subset of the grid
lo/hi = inclusive lo/hi bounds for brick of global grid cells I own
lo grid index = first grid pt >= fraclo bound
hi grid index = last grid pt < frachi bound
@ -420,12 +483,19 @@ void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
------------------------------------------------------------------------- */
void Grid3d::partition_grid(int ngrid, double fraclo, double frachi,
double shift, int &lo, int &hi)
double shift, int extra, int &lo, int &hi)
{
lo = static_cast<int> (fraclo * ngrid);
while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
hi = static_cast<int> (frachi * ngrid);
while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
if (extra == 0) {
lo = static_cast<int> (fraclo * ngrid);
while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
hi = static_cast<int> (frachi * ngrid);
while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
} else {
lo = static_cast<int> (fraclo * ngrid/zfactor);
while (1.0*lo + shift/ngrid < fraclo*ngrid) lo++;
hi = static_cast<int> (frachi * ngrid/zfactor);
while (1.0*hi + shift/ngrid >= frachi*ngrid) hi--;
}
}
/* ----------------------------------------------------------------------
@ -468,6 +538,8 @@ void Grid3d::ghost_grid()
double dxinv = nx / prd[0];
double dyinv = ny / prd[1];
double dzinv = nz / prd[2];
double dzinv_slab = nz / (prd[2] * zfactor);
int lo, hi;
lo = static_cast<int>((sublo[0]-dist[0]-boxlo[0]) * dxinv + shift_atom_lo + OFFSET) - OFFSET;
@ -486,10 +558,33 @@ void Grid3d::ghost_grid()
outzlo = MIN(lo - stencil_atom_lo, inzlo - stencil_grid_lo);
outzhi = MAX(hi + stencil_atom_hi, inzhi + stencil_grid_hi);
} else {
lo = static_cast<int>((sublo[2]-dist[2]-boxlo[2]) * dzinv_slab + shift_atom_lo + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[2]+dist[2]-boxlo[2]) * dzinv_slab + shift_atom_hi + OFFSET) - OFFSET;
outzlo = MIN(lo - stencil_atom_lo, inzlo - stencil_grid_lo);
outzhi = MAX(hi + stencil_atom_hi, inzhi + stencil_grid_hi);
}
// limit out xyz lo/hi indices to global grid for non-periodic dims
// if zextra = 1:
// adjust grid boundaries for processors at +z end,
// to include added empty grid cells between periodically repeating slabs
// in this case:
// want grid data forward communicated from +z proc to -z proc, but not vice versa
// want grid data reverse communicated from -z proc to +z proc, but not vice versa
// this is accomplished by inzhi = outzhi on +z end (no ghost cells)
// also insure no other procs use ghost cells beyond +z limit
if (zextra) {
if (layout != Comm::LAYOUT_TILED) {
if (comm->myloc[2] == comm->procgrid[2]-1) inzhi = outzhi = nz - 1;
} else {
if (comm->mysplit[2][1] == 1.0) inzhi = outzhi = nz - 1;
}
outzhi = MIN(outzhi,nz-1);
}
// limit out xyz lo/hi indices to global grid for non-periodic dims
// if zextra = 1 (e.g. PPPM), grid is still fully periodic
int *periodicity = domain->periodicity;
if (!periodicity[0]) {
@ -506,45 +601,6 @@ void Grid3d::ghost_grid()
outzlo = MAX(0,outzlo);
outzhi = MIN(nz-1,outzhi);
}
// NOTE: zperiod effects, also in owned grid cells ?
/*
// extent of zprd when 2d slab mode is selected
double zprd_slab = zprd*slab_volfactor;
// for slab PPPM, assign z grid as if it were not extended
nlo = static_cast<int> ((sublo[2]-dist[2]-boxlo[2]) *
nz_pppm/zprd_slab + shift) - OFFSET;
nhi = static_cast<int> ((subhi[2]+dist[2]-boxlo[2]) *
nz_pppm/zprd_slab + shift) - OFFSET;
nzlo_out = nlo + nlower;
nzhi_out = nhi + nupper;
// for slab PPPM, change the grid boundary for processors at +z end
// to include the empty volume between periodically repeating slabs
// for slab PPPM, want charge data communicated from -z proc to +z proc,
// but not vice versa, also want field data communicated from +z proc to
// -z proc, but not vice versa
// this is accomplished by nzhi_in = nzhi_out on +z end (no ghost cells)
// also insure no other procs use ghost cells beyond +z limit
// differnet logic for non-tiled vs tiled decomposition
if (slabflag == 1) {
if (comm->layout != Comm::LAYOUT_TILED) {
if (comm->myloc[2] == comm->procgrid[2]-1) nzhi_in = nzhi_out = nz_pppm - 1;
} else {
if (comm->mysplit[2][1] == 1.0) nzhi_in = nzhi_out = nz_pppm - 1;
}
nzhi_out = MIN(nzhi_out,nz_pppm-1);
}
*/
}
/* ----------------------------------------------------------------------
@ -1770,9 +1826,9 @@ int Grid3d::compute_overlap(int ghostflag, int *box, int *pbc, Overlap *&overlap
for (int k = kproclo; k <= kprochi; k++)
for (int j = jproclo; j <= jprochi; j++)
for (int i = iproclo; i <= iprochi; i++) {
partition_grid(nx,xsplit[i],xsplit[i+1],shift_grid,obox[0],obox[1]);
partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,obox[2],obox[3]);
partition_grid(nz,zsplit[k],zsplit[k+1],shift_grid,obox[4],obox[5]);
partition_grid(nx,xsplit[i],xsplit[i+1],shift_grid,0,obox[0],obox[1]);
partition_grid(ny,ysplit[j],ysplit[j+1],shift_grid,0,obox[2],obox[3]);
partition_grid(nz,zsplit[k],zsplit[k+1],shift_grid,zextra,obox[4],obox[5]);
if (noverlap_list == maxoverlap_list) grow_overlap();
overlap_list[noverlap_list].proc = grid2proc[i][j][k];

12
src/grid3d.h
View File

@ -23,6 +23,8 @@ class Grid3d : protected Pointers {
enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes
Grid3d(class LAMMPS *, MPI_Comm, int, int, int);
Grid3d(class LAMMPS *, MPI_Comm, int, int, int,
int, int, int, int, int, int, int, int, int, int, int, int);
~Grid3d();
void set_distance(double);
@ -31,7 +33,7 @@ class Grid3d : protected Pointers {
void set_shift_grid(double);
void set_shift_atom(double, double);
void set_zfactor(double);
void set_larger_grid(int, int, int, int, int, int);
void set_caller_grid(int, int, int, int, int, int);
void set_proc_neighs(int, int, int, int, int, int);
int identical(Grid3d *);
@ -59,7 +61,7 @@ class Grid3d : protected Pointers {
MPI_Comm gridcomm; // communicator for this class
// usually world, but MSM calls with subset
// inputs from caller
// input from caller
int nx, ny, nz; // size of global grid in all 3 dims
double maxdist; // distance owned atoms can move outside subdomain
@ -230,13 +232,11 @@ class Grid3d : protected Pointers {
// internal methods
// -------------------------------------------
void partition_grid(int, double, double, double, int &, int &);
void initialize();
void partition_grid(int, double, double, double, int, int &, int &);
void ghost_grid();
void extract_comm_info();
void store(int, int, int, int, int, int, int, int, int, int, int, int,
int, int, int, int, int, int, int, int, int, int, int, int);
void setup_comm_brick(int &, int &);
void setup_comm_tiled(int &, int &);
int ghost_adjacent_brick();