Merge pull request #2885 from lammps/grid-adjust

Grid adjust
This commit is contained in:
Axel Kohlmeyer 2021-08-23 21:39:35 -04:00 committed by GitHub
commit 43261c3a4f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 142 additions and 80 deletions

View File

@ -643,6 +643,7 @@ void MSM::allocate()
if (active_flag[n]) {
delete gc[n];
int **procneigh = procneigh_levels[n];
gc[n] = new GridComm(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],
@ -1194,20 +1195,13 @@ void MSM::set_grid_local()
for (int n=0; n<levels; n++) {
// global indices of MSM grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global MSM grid that I own without ghost cells
// 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
nxlo_in[n] = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx_msm[n]);
nxhi_in[n] = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx_msm[n]) - 1;
nylo_in[n] = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny_msm[n]);
nyhi_in[n] = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny_msm[n]) - 1;
nzlo_in[n] = static_cast<int> (comm->zsplit[comm->myloc[2]] * nz_msm[n]);
nzhi_in[n] = static_cast<int> (comm->zsplit[comm->myloc[2]+1] * nz_msm[n]) - 1;
// nlower,nupper = stencil size for mapping (interpolating) particles to MSM grid
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]);
nlower = -(order-1)/2;
nupper = order/2;
@ -1287,49 +1281,38 @@ void MSM::set_grid_local()
nzhi_out[n] = nhi + MAX(order,nzhi_direct);
// add extra grid points for non-periodic boundary conditions
// skip reset of lo/hi for procs who do not own any grid cells
if (domain->nonperiodic) {
if (!domain->xperiodic) {
if (nxlo_in[n] == 0)
nxlo_in[n] = alpha[n];
if (!domain->xperiodic && nxlo_in[n] <= nxhi_in[n]) {
if (nxlo_in[n] == 0) nxlo_in[n] = alpha[n];
nxlo_out[n] = MAX(nxlo_out[n],alpha[n]);
if (n == 0) nxlo_out_all = MAX(nxlo_out_all,alpha[0]);
if (nxhi_in[n] == nx_msm[n] - 1)
nxhi_in[n] = betax[n];
if (nxhi_in[n] == nx_msm[n] - 1) nxhi_in[n] = betax[n];
nxhi_out[n] = MIN(nxhi_out[n],betax[n]);
if (n == 0) nxhi_out_all = MIN(nxhi_out_all,betax[0]);
if (nxhi_in[n] < 0)
nxhi_in[n] = alpha[n] - 1;
}
if (!domain->yperiodic) {
if (nylo_in[n] == 0)
nylo_in[n] = alpha[n];
if (!domain->yperiodic && nylo_in[n] <= nyhi_in[n]) {
if (nylo_in[n] == 0) nylo_in[n] = alpha[n];
nylo_out[n] = MAX(nylo_out[n],alpha[n]);
if (n == 0) nylo_out_all = MAX(nylo_out_all,alpha[0]);
if (nyhi_in[n] == ny_msm[n] - 1)
nyhi_in[n] = betay[n];
if (nyhi_in[n] == ny_msm[n] - 1) nyhi_in[n] = betay[n];
nyhi_out[n] = MIN(nyhi_out[n],betay[n]);
if (n == 0) nyhi_out_all = MIN(nyhi_out_all,betay[0]);
if (nyhi_in[n] < 0)
nyhi_in[n] = alpha[n] - 1;
}
if (!domain->zperiodic) {
if (nzlo_in[n] == 0)
nzlo_in[n] = alpha[n];
if (!domain->zperiodic && nzlo_in[n] <= nzhi_in[n]) {
if (nzlo_in[n] == 0) nzlo_in[n] = alpha[n];
nzlo_out[n] = MAX(nzlo_out[n],alpha[n]);
if (n == 0) nzlo_out_all = MAX(nzlo_out_all,alpha[0]);
if (nzhi_in[n] == nz_msm[n] - 1)
nzhi_in[n] = betaz[n];
if (nzhi_in[n] == nz_msm[n] - 1) nzhi_in[n] = betaz[n];
nzhi_out[n] = MIN(nzhi_out[n],betaz[n]);
if (n == 0) nzhi_out_all = MIN(nzhi_out_all,betaz[0]);
if (nzhi_in[n] < 0)
nzhi_in[n] = alpha[n] - 1;
}
}

View File

@ -1320,34 +1320,12 @@ double PPPM::final_accuracy()
void PPPM::set_grid_local()
{
// global indices of PPPM grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global PPPM grid that I own without ghost cells
// for slab PPPM, assign z grid as if it were not extended
// both non-tiled and tiled proc layouts use 0-1 fractional sumdomain info
// 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
if (comm->layout != Comm::LAYOUT_TILED) {
nxlo_in = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx_pppm);
nxhi_in = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx_pppm) - 1;
nylo_in = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny_pppm);
nyhi_in = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny_pppm) - 1;
nzlo_in = static_cast<int>
(comm->zsplit[comm->myloc[2]] * nz_pppm/slab_volfactor);
nzhi_in = static_cast<int>
(comm->zsplit[comm->myloc[2]+1] * nz_pppm/slab_volfactor) - 1;
} else {
nxlo_in = static_cast<int> (comm->mysplit[0][0] * nx_pppm);
nxhi_in = static_cast<int> (comm->mysplit[0][1] * nx_pppm) - 1;
nylo_in = static_cast<int> (comm->mysplit[1][0] * ny_pppm);
nyhi_in = static_cast<int> (comm->mysplit[1][1] * ny_pppm) - 1;
nzlo_in = static_cast<int> (comm->mysplit[2][0] * nz_pppm/slab_volfactor);
nzhi_in = static_cast<int> (comm->mysplit[2][1] * nz_pppm/slab_volfactor) - 1;
}
comm->partition_grid(nx_pppm,ny_pppm,nz_pppm,slab_volfactor,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
// nlower,nupper = stencil size for mapping particles to PPPM grid
@ -1367,7 +1345,7 @@ void PPPM::set_grid_local()
// effectively nlo_in,nhi_in + ghost cells
// nlo,nhi = index of global grid pt to "lower left" of smallest/largest
// position a particle in my box can be at
// dist[3] = particle position bound = subbox + skin/2.0 + qdist
// dist[3] = max particle position outside subbox = skin/2.0 + qdist
// qdist = offset due to TIP4P fictitious charge
// convert to triclinic if necessary
// nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping

View File

@ -2700,22 +2700,12 @@ void PPPMDisp::set_fft_parameters(int& nx_p, int& ny_p, int& nz_p,
int& ng, int& nf, int& nfb,
double& sft, double& sftone, int& ord)
{
// global indices of PPPM grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global PPPM grid that I own without ghost cells
// for slab PPPM, assign z grid as if it were not extended
nxlo_i = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx_p);
nxhi_i = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx_p) - 1;
nylo_i = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny_p);
nyhi_i = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny_p) - 1;
nzlo_i = static_cast<int>
(comm->zsplit[comm->myloc[2]] * nz_p/slab_volfactor);
nzhi_i = static_cast<int>
(comm->zsplit[comm->myloc[2]+1] * nz_p/slab_volfactor) - 1;
// partition global grid across procs
// n xyz lo/hi i = lower/upper bounds of global grid this proc owns
// indices range from 0 to N-1 inclusive in each dim
comm->partition_grid(nx_p,ny_p,nz_p,slab_volfactor,
nxlo_i,nxhi_i,nylo_i,nyhi_i,nzlo_i,nzhi_i);
// nlow,nupp = stencil size for mapping particles to PPPM grid

View File

@ -832,6 +832,103 @@ int Comm::binary(double value, int n, double *vec)
return index;
}
/* ----------------------------------------------------------------------
partition a global regular grid into one brick-shaped sub-grid per proc
if grid point is inside my sub-domain I own it,
this includes sub-domain lo boundary but excludes hi boundary
nx,ny,nz = extent of global grid
indices into the global grid range from 0 to N-1 in each dim
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
nxyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own
if proc owns no grid cells in a dim, then nlo > nhi
special case: 2 procs share boundary which a grid point is exactly on
2 equality if tests insure a consistent decision as to which proc owns it
------------------------------------------------------------------------- */
void Comm::partition_grid(int nx, int ny, int nz, double zfactor,
int &nxlo, int &nxhi, int &nylo, int &nyhi,
int &nzlo, int &nzhi)
{
double xfraclo,xfrachi,yfraclo,yfrachi,zfraclo,zfrachi;
if (layout != LAYOUT_TILED) {
xfraclo = xsplit[myloc[0]];
xfrachi = xsplit[myloc[0]+1];
yfraclo = ysplit[myloc[1]];
yfrachi = ysplit[myloc[1]+1];
zfraclo = zsplit[myloc[2]];
zfrachi = zsplit[myloc[2]+1];
} else {
xfraclo = mysplit[0][0];
xfrachi = mysplit[0][1];
yfraclo = mysplit[1][0];
yfrachi = mysplit[1][1];
zfraclo = mysplit[2][0];
zfrachi = mysplit[2][1];
}
nxlo = static_cast<int> (xfraclo * nx);
if (1.0*nxlo != xfraclo*nx) nxlo++;
nxhi = static_cast<int> (xfrachi * nx);
if (1.0*nxhi == xfrachi*nx) nxhi--;
nylo = static_cast<int> (yfraclo * ny);
if (1.0*nylo != yfraclo*ny) nylo++;
nyhi = static_cast<int> (yfrachi * ny);
if (1.0*nyhi == yfrachi*ny) nyhi--;
if (zfactor == 0.0) {
nzlo = static_cast<int> (zfraclo * nz);
if (1.0*nzlo != zfraclo*nz) nzlo++;
nzhi = static_cast<int> (zfrachi * nz);
if (1.0*nzhi == zfrachi*nz) nzhi--;
} else {
nzlo = static_cast<int> (zfraclo * nz/zfactor);
if (1.0*nzlo != zfraclo*nz) nzlo++;
nzhi = static_cast<int> (zfrachi * nz/zfactor);
if (1.0*nzhi == zfrachi*nz) nzhi--;
}
// OLD code
// could sometimes map grid points slightly outside a proc to the proc
/*
if (layout != LAYOUT_TILED) {
nxlo = static_cast<int> (xsplit[myloc[0]] * nx);
nxhi = static_cast<int> (xsplit[myloc[0]+1] * nx) - 1;
nylo = static_cast<int> (ysplit[myloc[1]] * ny);
nyhi = static_cast<int> (ysplit[myloc[1]+1] * ny) - 1;
if (zfactor == 0.0) {
nzlo = static_cast<int> (zsplit[myloc[2]] * nz);
nzhi = static_cast<int> (zsplit[myloc[2]+1] * nz) - 1;
} else {
nzlo = static_cast<int> (zsplit[myloc[2]] * nz/zfactor);
nzhi = static_cast<int> (zsplit[myloc[2]+1] * nz/zfactor) - 1;
}
} else {
nxlo = static_cast<int> (mysplit[0][0] * nx);
nxhi = static_cast<int> (mysplit[0][1] * nx) - 1;
nylo = static_cast<int> (mysplit[1][0] * ny);
nyhi = static_cast<int> (mysplit[1][1] * ny) - 1;
if (zfactor == 0.0) {
nzlo = static_cast<int> (mysplit[2][0] * nz);
nzhi = static_cast<int> (mysplit[2][1] * nz) - 1;
} else {
nzlo = static_cast<int> (mysplit[2][0] * nz/zfactor);
nzhi = static_cast<int> (mysplit[2][1] * nz/zfactor) - 1;
}
}
*/
}
/* ----------------------------------------------------------------------
communicate inbuf around full ring of processors with messtag
nbytes = size of inbuf = n datums * nper bytes

View File

@ -105,6 +105,11 @@ class Comm : protected Pointers {
virtual void coord2proc_setup() {}
virtual int coord2proc(double *, int &, int &, int &);
// partition a global regular grid by proc sub-domains
void partition_grid(int, int, int, double,
int &, int &, int &, int &, int &, int &);
// memory usage
virtual double memory_usage() = 0;
@ -117,6 +122,7 @@ class Comm : protected Pointers {
int statflag = 0);
// extract data useful to other classes
virtual void *extract(const char *, int &) { return nullptr; }
protected:

View File

@ -42,38 +42,46 @@ class Domain : protected Pointers {
int tiltsmall; // 1 if limit tilt, else 0
// orthogonal box
double xprd, yprd, zprd; // global box dimensions
double xprd_half, yprd_half, zprd_half; // half dimensions
double prd[3]; // array form of dimensions
double prd_half[3]; // array form of half dimensions
// triclinic box
// xprd,xprd_half,prd,prd_half =
// same as if untilted
// xyzprd,xyzprd_half and prd,prd_half = same as if untilted
double prd_lamda[3]; // lamda box = (1,1,1)
double prd_half_lamda[3]; // lamda half box = (0.5,0.5,0.5)
double boxlo[3], boxhi[3]; // orthogonal box global bounds
// orthogonal box global bounds
double boxlo[3], boxhi[3];
// triclinic box
// boxlo/hi = same as if untilted
double boxlo_lamda[3], boxhi_lamda[3]; // lamda box = (0,1)
double boxlo_bound[3], boxhi_bound[3]; // bounding box of tilted domain
double corners[8][3]; // 8 corner points
// orthogonal box & triclinic box
double minxlo, minxhi; // minimum size of global box
double minylo, minyhi; // when shrink-wrapping
double minzlo, minzhi; // tri only possible for non-skew dims
// orthogonal box
double sublo[3], subhi[3]; // sub-box bounds on this proc
// triclinic box
// sublo/hi = undefined
double sublo_lamda[3], subhi_lamda[3]; // bounds of subbox in lamda
// triclinic box
double xy, xz, yz; // 3 tilt factors
double h[6], h_inv[6]; // shape matrix in Voigt ordering
// Voigt = xx,yy,zz,yz,xz,xy