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

more classes use Grid3d

Browse Source
This commit is contained in:
Steve Plimpton
2022-11-01 15:53:39 -06:00
parent 335bacebb7
commit 4c29457351
14 changed files with 440 additions and 393 deletions
Download Patch File Download Diff File Expand all files Collapse all files

114
src/AMOEBA/amoeba_convolution.cpp
View File

@ -104,82 +104,29 @@ void AmoebaConvolution::reset_grid()
void AmoebaConvolution::allocate_grid()
{
// global indices of grid range from 0 to N-1
// nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of
// global grid that I own without ghost cells
// both non-tiled and tiled proc layouts use 0-1 fractional subdomain info
if (comm->layout != Comm::LAYOUT_TILED) {
nxlo_in = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx);
nxhi_in = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx) - 1;
nylo_in = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny);
nyhi_in = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny) - 1;
nzlo_in = static_cast<int> (comm->zsplit[comm->myloc[2]] * nz);
nzhi_in = static_cast<int> (comm->zsplit[comm->myloc[2]+1] * nz) - 1;
} else {
nxlo_in = static_cast<int> (comm->mysplit[0][0] * nx);
nxhi_in = static_cast<int> (comm->mysplit[0][1] * nx) - 1;
nylo_in = static_cast<int> (comm->mysplit[1][0] * ny);
nyhi_in = static_cast<int> (comm->mysplit[1][1] * ny) - 1;
nzlo_in = static_cast<int> (comm->mysplit[2][0] * nz);
nzhi_in = static_cast<int> (comm->mysplit[2][1] * nz) - 1;
}
// maxdist = max distance outside of subbox my owned atom may be
// nlower,nupper = stencil size for mapping particles to FFT grid
double maxdist = 0.5*neighbor->skin;
int nlower = -(order-1)/2;
int nupper = order/2;
// nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of
// global grid that my particles can contribute charge to
// effectively nlo_in,nhi_in + ghost cells
// nlo,nhi = global coords of 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
// convert to triclinic if necessary
// nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping
// Grid3d determines my owned + ghost grid cells
// ghost cell extent depends on maxdist, nlower, nupper
gc = new Grid3d(lmp,world,nx,ny,nz);
gc->set_distance(maxdist);
gc->set_stencil_atom(-nlower,nupper);
double *prd,*boxlo,*sublo,*subhi;
int triclinic = domain->triclinic;
gc->setup_grid(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
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;
}
int nqty = flag3d ? 1 : 2;
int ngc_buf1,ngc_buf2;
double xprd = prd[0];
double yprd = prd[1];
double zprd = prd[2];
double dist[3] = {0.0,0.0,0.0};
double cuthalf = 0.5*neighbor->skin;
if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf;
else MathExtra::tribbox(domain->h,cuthalf,&dist[0]);
int nlo,nhi;
nlo = static_cast<int> ((sublo[0]-dist[0]-boxlo[0]) * nx/xprd);
nhi = static_cast<int> ((subhi[0]+dist[0]-boxlo[0]) * nx/xprd);
nxlo_out = nlo + nlower;
nxhi_out = nhi + nupper;
nlo = static_cast<int> ((sublo[1]-dist[1]-boxlo[1]) * ny/yprd);
nhi = static_cast<int> ((subhi[1]+dist[1]-boxlo[1]) * ny/yprd);
nylo_out = nlo + nlower;
nyhi_out = nhi + nupper;
nlo = static_cast<int> ((sublo[2]-dist[2]-boxlo[2]) * nz/zprd);
nhi = static_cast<int> ((subhi[2]+dist[2]-boxlo[2]) * nz/zprd);
nzlo_out = nlo + nlower;
nzhi_out = nhi + nupper;
gc->setup_comm(ngc_buf1,ngc_buf2);
memory->create(gc_buf1,nqty*ngc_buf1,"amoeba:gc_buf1");
memory->create(gc_buf2,nqty*ngc_buf2,"amoeba:gc_buf2");
// x-pencil decomposition of FFT mesh
// global indices range from 0 to N-1
@ -210,7 +157,7 @@ void AmoebaConvolution::allocate_grid()
nzlo_fft = me_z*nz/npez_fft;
nzhi_fft = (me_z+1)*nz/npez_fft - 1;
// grid sizes
// FFT grid sizes
// nbrick_owned = owned grid points in brick decomp
// nbrick_ghosts = owned + ghost grid points in grid decomp
// nfft_owned = owned grid points in FFT decomp
@ -226,7 +173,7 @@ void AmoebaConvolution::allocate_grid()
ngrid_either = MAX(nbrick_owned,nfft_owned);
// instantiate FFT, Grid3d, and Remap
// instantiate FFT and Remap
int tmp;
@ -240,11 +187,6 @@ void AmoebaConvolution::allocate_grid()
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
0,0,&tmp,0);
gc = new Grid3d(lmp,world,nx,ny,nz,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
int nqty = flag3d ? 1 : 2;
remap = new Remap(lmp,world,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
@ -266,12 +208,6 @@ void AmoebaConvolution::allocate_grid()
memory->create(grid_fft,ngrid_either,"amoeba:grid_fft");
memory->create(cfft,2*ngrid_either,"amoeba:cfft");
int ngc_buf1,ngc_buf2;
gc->setup(ngc_buf1,ngc_buf2);
memory->create(gc_buf1,nqty*ngc_buf1,"amoeba:gc_buf1");
memory->create(gc_buf2,nqty*ngc_buf2,"amoeba:gc_buf2");
memory->create(remap_buf,nqty*nfft_owned,"amoeba:remap_buf");
}
@ -279,18 +215,20 @@ void AmoebaConvolution::allocate_grid()
void AmoebaConvolution::deallocate_grid()
{
delete gc;
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
delete fft1;
delete fft2;
delete remap;
memory->destroy3d_offset(grid_brick,nzlo_out,nylo_out,nxlo_out);
memory->destroy4d_offset_last(cgrid_brick,nzlo_out,nylo_out,nxlo_out);
memory->destroy(grid_fft);
memory->destroy(cfft);
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
memory->destroy(remap_buf);
delete fft1;
delete fft2;
delete gc;
delete remap;
}
/* ----------------------------------------------------------------------

4
src/AMOEBA/amoeba_kspace.cpp
View File

@ -217,6 +217,10 @@ void PairAmoeba::bspline_fill()
// NOTE: could subtract off nlpts to start with
// NOTE: this is place to check that stencil size does not
// go out of bounds relative to igrid for a proc's sub-domain
// similar to PPPM::particle_map()
// subtracting eps is strange, and could mess up the check
// better to check here than in methods like grid_mpole()
// but check needs to be valid for all KSpace terms
// NOTE: could convert x -> lamda for entire set of Nlocal atoms
domain->x2lamda(x[i],lamda);

17
src/EXTRA-FIX/fix_ttm_grid.cpp
View File

@ -485,10 +485,11 @@ void FixTTMGrid::reset_grid()
int tmp[12];
double maxdist = 0.5 * neighbor->skin;
Grid3d *gridnew = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid,
maxdist, 1, 0.5,
tmp[0],tmp[1],tmp[2],tmp[3],tmp[4],tmp[5],
tmp[6],tmp[7],tmp[8],tmp[9],tmp[10],tmp[11]);
Grid3d *gridnew = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
gridnew->set_distance(maxdist);
gridnew->set_stencil_grid(1,1);
gridnew->setup_grid(tmp[0],tmp[1],tmp[2],tmp[3],tmp[4],tmp[5],
tmp[6],tmp[7],tmp[8],tmp[9],tmp[10],tmp[11]);
if (grid->identical(gridnew)) {
delete gridnew;
@ -632,9 +633,11 @@ void FixTTMGrid::allocate_grid()
{
double maxdist = 0.5 * neighbor->skin;
grid = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, maxdist, 1, 0.5,
nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
grid = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
grid->set_distance(maxdist);
grid->set_stencil_grid(1,1);
grid->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
ngridown = (nxhi_in - nxlo_in + 1) * (nyhi_in - nylo_in + 1) *
(nzhi_in - nzlo_in + 1);

4
src/KSPACE/msm.cpp
View File

@ -604,7 +604,7 @@ void MSM::allocate()
nxlo_out[0],nxhi_out[0],nylo_out[0],
nyhi_out[0],nzlo_out[0],nzhi_out[0]);
gcall->setup(ngcall_buf1,ngcall_buf2);
gcall->setup_comm(ngcall_buf1,ngcall_buf2);
npergrid = 1;
memory->destroy(gcall_buf1);
memory->destroy(gcall_buf2);
@ -636,7 +636,7 @@ void MSM::allocate()
procneigh[0][0],procneigh[0][1],procneigh[1][0],
procneigh[1][1],procneigh[2][0],procneigh[2][1]);
gc[n]->setup(ngc_buf1[n],ngc_buf2[n]);
gc[n]->setup_comm(ngc_buf1[n],ngc_buf2[n]);
npergrid = 1;
memory->destroy(gc_buf1[n]);
memory->destroy(gc_buf2[n]);

166
src/KSPACE/pppm.cpp
View File

@ -292,7 +292,7 @@ void PPPM::init()
// or overlap is allowed, then done
// else reduce order and try again
Grid3d *gctmp = nullptr;
gc = nullptr;
int iteration = 0;
while (order >= minorder) {
@ -305,23 +305,28 @@ void PPPM::init()
set_grid_local();
if (overlap_allowed) break;
gctmp = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm);
gc->set_distance(0.5*neighbor->skin + qdist);
gc->set_stencil_atom(-nlower,nupper);
gc->set_shift_atom(shiftatom);
gc->set_zfactor(slab_volfactor);
gc->setup_grid(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
int tmp1,tmp2;
gctmp->setup(tmp1,tmp2);
if (gctmp->ghost_adjacent()) break;
delete gctmp;
gc->setup_comm(tmp1,tmp2);
if (gc->ghost_adjacent()) break;
delete gc;
order--;
iteration++;
}
if (order < minorder) error->all(FLERR,"PPPM order < minimum allowed order");
if (!overlap_allowed && !gctmp->ghost_adjacent())
if (!overlap_allowed && !gc->ghost_adjacent())
error->all(FLERR,"PPPM grid stencil extends beyond nearest neighbor processor");
if (gctmp) delete gctmp;
if (gc) delete gc;
// adjust g_ewald
@ -739,6 +744,29 @@ void PPPM::compute(int eflag, int vflag)
void PPPM::allocate()
{
// create ghost grid object for rho and electric field communication
// returns local owned and ghost grid bounds
// setup communication patterns and buffers
gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm);
gc->set_distance(0.5*neighbor->skin + qdist);
gc->set_stencil_atom(-nlower,nupper);
gc->set_shift_atom(shiftatom);
gc->set_zfactor(slab_volfactor);
gc->setup_grid(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc->setup_comm(ngc_buf1,ngc_buf2);
if (differentiation_flag) npergrid = 1;
else npergrid = 3;
memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1");
memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2");
// allocate distributed grid data
memory->create3d_offset(density_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"pppm:density_brick");
@ -809,21 +837,6 @@ void PPPM::allocate()
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
1,0,0,FFT_PRECISION,collective_flag);
// create ghost grid object for rho and electric field communication
// also create 2 bufs for ghost grid cell comm, passed to Grid3d methods
gc = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc->setup(ngc_buf1,ngc_buf2);
if (differentiation_flag) npergrid = 1;
else npergrid = 3;
memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1");
memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2");
}
/* ----------------------------------------------------------------------
@ -832,6 +845,10 @@ void PPPM::allocate()
void PPPM::deallocate()
{
delete gc;
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
memory->destroy3d_offset(density_brick,nzlo_out,nylo_out,nxlo_out);
if (differentiation_flag == 1) {
@ -874,9 +891,6 @@ void PPPM::deallocate()
delete fft1;
delete fft2;
delete remap;
delete gc;
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
}
/* ----------------------------------------------------------------------
@ -1309,85 +1323,41 @@ double PPPM::final_accuracy()
}
/* ----------------------------------------------------------------------
set local subset of PPPM/FFT grid that I own
n xyz lo/hi in = 3d brick that I own (inclusive)
n xyz lo/hi out = 3d brick + ghost cells in 6 directions (inclusive)
n xyz lo/hi fft = FFT columns that I own (all of x dim, 2d decomp in yz)
set params which determine which owned and ghost cells this proc owns
Grid3d will use these params to partition grid
also partition FFT grid
n xyz lo/hi fft = FFT columns that I own (all of x dim, 2d decomp in yz)
------------------------------------------------------------------------- */
void PPPM::set_grid_local()
{
// 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
// shift values for particle <-> grid mapping
// add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
comm->partition_grid(nx_pppm,ny_pppm,nz_pppm,slab_volfactor,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
if (order % 2) shiftatom = 0.5;
else shiftatom = 0.0;
shift = OFFSET + shiftatom;
if (order % 2) shiftone = 0.0;
else shiftone = 0.5;
// nlower,nupper = stencil size for mapping particles to PPPM grid
nlower = -(order-1)/2;
nupper = order/2;
// shift values for particle <-> grid mapping
// add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
if (order % 2) shift = OFFSET + 0.5;
else shift = OFFSET;
if (order % 2) shiftone = 0.0;
else shiftone = 0.5;
// nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of
// global PPPM grid that my particles can contribute charge to
// 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] = 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
// for slab PPPM, assign z grid as if it were not extended
// NOTE: still to do: stagger and zperiod effects
double *prd,*sublo,*subhi;
/*
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];
double yprd = prd[1];
double zprd = prd[2];
// extent of zprd when 2d slab mode is selected
double zprd_slab = zprd*slab_volfactor;
double dist[3] = {0.0,0.0,0.0};
double cuthalf = 0.5*neighbor->skin + qdist;
if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf;
else MathExtra::tribbox(domain->h,cuthalf,&dist[0]);
int nlo,nhi;
nlo = nhi = 0;
nlo = static_cast<int> ((sublo[0]-dist[0]-boxlo[0]) *
nx_pppm/xprd + shift) - OFFSET;
nhi = static_cast<int> ((subhi[0]+dist[0]-boxlo[0]) *
nx_pppm/xprd + shift) - OFFSET;
nxlo_out = nlo + nlower;
nxhi_out = nhi + nupper;
nlo = static_cast<int> ((sublo[1]-dist[1]-boxlo[1]) *
ny_pppm/yprd + shift) - OFFSET;
nhi = static_cast<int> ((subhi[1]+dist[1]-boxlo[1]) *
ny_pppm/yprd + shift) - OFFSET;
nylo_out = nlo + nlower;
nyhi_out = nhi + nupper;
// 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;
@ -1396,6 +1366,7 @@ void PPPM::set_grid_local()
nzlo_out = nlo + nlower;
nzhi_out = nhi + nupper;
if (stagger_flag) {
nxhi_out++;
nyhi_out++;
@ -1420,6 +1391,11 @@ void PPPM::set_grid_local()
nzhi_out = MIN(nzhi_out,nz_pppm-1);
}
*/
// x-pencil decomposition of FFT mesh
// global indices range from 0 to N-1
// each proc owns entire x-dimension, clumps of columns in y,z dimensions
@ -1446,18 +1422,22 @@ void PPPM::set_grid_local()
nzlo_fft = me_z*nz_pppm/npez_fft;
nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1;
// nfft = FFT points in x-pencil FFT decomposition on this proc
nfft = (nxhi_fft-nxlo_fft+1) * (nyhi_fft-nylo_fft+1) *
(nzhi_fft-nzlo_fft+1);
// ngrid = count of PPPM grid pts owned by this proc, including ghosts
ngrid = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) *
(nzhi_out-nzlo_out+1);
// count of FFT grids pts owned by this proc, without ghosts
// nfft = FFT points in x-pencil FFT decomposition on this proc
// nfft_brick = FFT points in 3d brick-decomposition on this proc
// nfft_both = greater of 2 values
nfft = (nxhi_fft-nxlo_fft+1) * (nyhi_fft-nylo_fft+1) *
(nzhi_fft-nzlo_fft+1);
int nfft_brick = (nxhi_in-nxlo_in+1) * (nyhi_in-nylo_in+1) *
(nzhi_in-nzlo_in+1);
nfft_both = MAX(nfft,nfft_brick);
@ -1908,7 +1888,7 @@ void PPPM::make_rho()
// loop over my charges, add their contribution to nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// (mx,my,mz) = global indices of moving stencil pt
double *q = atom->q;
double **x = atom->x;

2
src/KSPACE/pppm.h
View File

@ -48,7 +48,7 @@ class PPPM : public KSpace {
double volume;
double delxinv, delyinv, delzinv, delvolinv;
double h_x, h_y, h_z;
double shift, shiftone;
double shift, shiftone, shiftatom;
int peratom_allocate_flag;
int nxlo_in, nylo_in, nzlo_in, nxhi_in, nyhi_in, nzhi_in;

70
src/KSPACE/pppm_dipole.cpp
View File

@ -28,6 +28,7 @@
#include "math_const.h"
#include "math_special.h"
#include "memory.h"
#include "neighbor.h"
#include "pair.h"
#include "remap_wrap.h"
#include "update.h"
@ -188,7 +189,7 @@ void PPPMDipole::init()
// or overlap is allowed, then done
// else reduce order and try again
Grid3d *gctmp = nullptr;
gc_dipole = nullptr;
int iteration = 0;
while (order >= minorder) {
@ -201,23 +202,28 @@ void PPPMDipole::init()
set_grid_local();
if (overlap_allowed) break;
gctmp = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm);
gc_dipole->set_distance(0.5*neighbor->skin + qdist);
gc_dipole->set_stencil_atom(-nlower,nupper);
gc_dipole->set_shift_atom(shiftatom);
gc_dipole->set_zfactor(slab_volfactor);
gc_dipole->setup_grid(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
int tmp1,tmp2;
gctmp->setup(tmp1,tmp2);
if (gctmp->ghost_adjacent()) break;
delete gctmp;
gc_dipole->setup_comm(tmp1,tmp2);
if (gc_dipole->ghost_adjacent()) break;
delete gc_dipole;
order--;
iteration++;
}
if (order < minorder) error->all(FLERR,"PPPMDipole order < minimum allowed order");
if (!overlap_allowed && !gctmp->ghost_adjacent())
if (!overlap_allowed && !gc_dipole->ghost_adjacent())
error->all(FLERR,"PPPMDipole grid stencil extends beyond nearest neighbor processor");
if (gctmp) delete gctmp;
if (gc_dipole) delete gc_dipole;
// adjust g_ewald
@ -436,7 +442,7 @@ void PPPMDipole::compute(int eflag, int vflag)
particle_map();
make_rho_dipole();
// all procs communicate density values from their ghost cells
// to fully sum contribution in their 3d bricks
// remap from 3d decomposition to FFT decomposition
@ -528,6 +534,28 @@ void PPPMDipole::compute(int eflag, int vflag)
void PPPMDipole::allocate()
{
// create ghost grid object for rho and electric field communication
// returns local owned and ghost grid bounds
// setup communication patterns and buffers
gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm);
gc_dipole->set_distance(0.5*neighbor->skin + qdist);
gc_dipole->set_stencil_atom(-nlower,nupper);
gc_dipole->set_shift_atom(shiftatom);
gc_dipole->set_zfactor(slab_volfactor);
gc_dipole->setup_grid(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc_dipole->setup_comm(ngc_buf1,ngc_buf2);
npergrid = 9;
memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1");
memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2");
// allocate distributed grid data
memory->create3d_offset(densityx_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
nxlo_out,nxhi_out,"pppm_dipole:densityx_brick_dipole");
memory->create3d_offset(densityy_brick_dipole,nzlo_out,nzhi_out,nylo_out,nyhi_out,
@ -601,20 +629,6 @@ void PPPMDipole::allocate()
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft,
1,0,0,FFT_PRECISION,collective_flag);
// create ghost grid object for rho and electric field communication
// also create 2 bufs for ghost grid cell comm, passed to Grid3d methods
gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc_dipole->setup(ngc_buf1,ngc_buf2);
npergrid = 9;
memory->create(gc_buf1,npergrid*ngc_buf1,"pppm:gc_buf1");
memory->create(gc_buf2,npergrid*ngc_buf2,"pppm:gc_buf2");
}
/* ----------------------------------------------------------------------
@ -623,6 +637,10 @@ void PPPMDipole::allocate()
void PPPMDipole::deallocate()
{
delete gc_dipole;
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
memory->destroy3d_offset(densityx_brick_dipole,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(densityy_brick_dipole,nzlo_out,nylo_out,nxlo_out);
memory->destroy3d_offset(densityz_brick_dipole,nzlo_out,nylo_out,nxlo_out);
@ -645,7 +663,9 @@ void PPPMDipole::deallocate()
memory->destroy(work3);
memory->destroy(work4);
delete gc_dipole;
delete fft1;
delete fft2;
delete remap;
}
/* ----------------------------------------------------------------------

24
src/KSPACE/pppm_dipole_spin.cpp
View File

@ -26,6 +26,7 @@
#include "grid3d.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "pair.h"
#include "update.h"
@ -173,7 +174,7 @@ void PPPMDipoleSpin::init()
// or overlap is allowed, then done
// else reduce order and try again
Grid3d *gctmp = nullptr;
gc_dipole = nullptr;
int iteration = 0;
while (order >= minorder) {
@ -186,23 +187,28 @@ void PPPMDipoleSpin::init()
set_grid_local();
if (overlap_allowed) break;
gctmp = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm,
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc_dipole = new Grid3d(lmp,world,nx_pppm,ny_pppm,nz_pppm);
gc_dipole->set_distance(0.5*neighbor->skin + qdist);
gc_dipole->set_stencil_atom(-nlower,nupper);
gc_dipole->set_shift_atom(shiftatom);
gc_dipole->set_zfactor(slab_volfactor);
gc_dipole->setup_grid(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
int tmp1,tmp2;
gctmp->setup(tmp1,tmp2);
if (gctmp->ghost_adjacent()) break;
delete gctmp;
gc_dipole->setup_comm(tmp1,tmp2);
if (gc_dipole->ghost_adjacent()) break;
delete gc_dipole;
order--;
iteration++;
}
if (order < minorder) error->all(FLERR,"PPPMDipoleSpin order < minimum allowed order");
if (!overlap_allowed && !gctmp->ghost_adjacent())
if (!overlap_allowed && !gc_dipole->ghost_adjacent())
error->all(FLERR,"PPPMDipoleSpin grid stencil extends beyond nearest neighbor processor");
if (gctmp) delete gctmp;
if (gc_dipole) delete gc_dipole;
// adjust g_ewald

12
src/KSPACE/pppm_disp.cpp
View File

@ -446,7 +446,7 @@ void PPPMDisp::init()
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
int tmp1,tmp2;
gctmp->setup(tmp1,tmp2);
gctmp->setup_comm(tmp1,tmp2);
if (gctmp->ghost_adjacent()) break;
delete gctmp;
@ -526,7 +526,7 @@ void PPPMDisp::init()
nzlo_out_6,nzhi_out_6);
int tmp1,tmp2;
gctmp->setup(tmp1,tmp2);
gctmp->setup_comm(tmp1,tmp2);
if (gctmp->ghost_adjacent()) break;
delete gctmp;
@ -1754,7 +1754,7 @@ void _noopt PPPMDisp::allocate()
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out);
gc->setup(ngc_buf1,ngc_buf2);
gc->setup_comm(ngc_buf1,ngc_buf2);
if (differentiation_flag) npergrid = 1;
else npergrid = 3;
@ -1838,7 +1838,7 @@ void _noopt PPPMDisp::allocate()
nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6,
nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6);
gc6->setup(ngc6_buf1,ngc6_buf2);
gc6->setup_comm(ngc6_buf1,ngc6_buf2);
if (differentiation_flag) npergrid6 = 1;
else npergrid6 = 3;
@ -2001,7 +2001,7 @@ void _noopt PPPMDisp::allocate()
nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6,
nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6);
gc6->setup(ngc6_buf1,ngc6_buf2);
gc6->setup_comm(ngc6_buf1,ngc6_buf2);
if (differentiation_flag) npergrid6 = 7;
else npergrid6 = 21;
@ -2088,7 +2088,7 @@ void _noopt PPPMDisp::allocate()
nxlo_in_6,nxhi_in_6,nylo_in_6,nyhi_in_6,nzlo_in_6,nzhi_in_6,
nxlo_out_6,nxhi_out_6,nylo_out_6,nyhi_out_6,nzlo_out_6,nzhi_out_6);
gc6->setup(ngc6_buf1,ngc6_buf2);
gc6->setup_comm(ngc6_buf1,ngc6_buf2);
if (differentiation_flag) npergrid6 = 1*nsplit_alloc;
else npergrid6 = 3*nsplit_alloc;

6
src/compute_property_grid.cpp
View File

@ -252,9 +252,9 @@ void ComputePropertyGrid::allocate_grid()
ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1);
} else {
grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, 0.0, 0, 0.0, nxlo_in, nxhi_in, nylo_in,
nyhi_in, nzlo_in, nzhi_in, nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out,
nzhi_out);
grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
grid3d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
if (nvalues == 1)
memory->create3d_offset(vec3d, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out, nxhi_out,
"property/grid:vec3d");

2
src/dump_grid.cpp
View File

@ -509,7 +509,7 @@ int DumpGrid::count()
grid3d = (Grid3d *) compute[field2index[0]]->get_grid_by_index(field2grid[0]);
else if (field2source[0] == FIX)
grid3d = (Grid3d *) fix[field2index[0]]->get_grid_by_index(field2grid[0]);
grid3d->get_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
grid3d->get_bounds_owned(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
}
// invoke Computes for per-grid quantities

8
src/fix_ave_grid.cpp
View File

@ -411,10 +411,10 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
ngridout = (nxhi_out - nxlo_out + 1) * (nyhi_out - nylo_out + 1);
} else {
grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, maxdist, 0, 0.5,
nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out,
nzlo_out, nzhi_out);
grid3d = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid);
grid3d->set_distance(maxdist);
grid3d->setup_grid(nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
// ngrid_buf12 converted to nvalues + count

367
src/grid3d.cpp
View File

@ -40,132 +40,29 @@ static constexpr int OFFSET = 16384;
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
constructor called by all classes except PPPM and MSM
comm_caller = caller's communicator
nx,ny,nz caller = size of global grid
maxdist = max distance outside of proc domain a particle will be
extra = additional ghost grid pts needed in each dim, e.g. for stencil
shift_caller = 0.0 for grid pt in lower-left corner of grid cell,
0.5 for center
return:
i xyz lohi = portion of global grid this proc owns, 0 <= index < N
o xyz lohi = owned + ghost grid cells needed in all directions
for periodic dims, o indices can be < 0 or >= N
for non-periodic dims, o indices will be >= 0 and < N
since no grid comm is done across non-periodic boundaries
NOTE: allow zfactor to be a calling arg for PPPM ?
constructor to create a 3d distributed grid
gcomm = caller's communicator
gnx,gny,gnz = global grid size
------------------------------------------------------------------------- */
Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm comm_caller,
int nx_caller, int ny_caller, int nz_caller,
double maxdist, int extra, double shift_caller,
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)
Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int gnx, int gny, int gnz) :
Pointers(lmp)
{
// store commnicator and global grid size
// set layout mode
gridcomm = comm_caller;
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = nx_caller;
ny = ny_caller;
nz = nz_caller;
nx = gnx;
ny = gny;
nz = gnz;
// default settings, can be overridden by set() methods
// NOTE: hardwire shift = 0.5 ?
shift = shift_caller;
// define owned grid cells for each proc
// extend bounds with ghost grid cells in each direction
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,inxlo,inxhi);
fraclo = comm->ysplit[comm->myloc[1]];
frachi = comm->ysplit[comm->myloc[1]+1];
partition_grid(ny,fraclo,frachi,shift,inylo,inyhi);
fraclo = comm->zsplit[comm->myloc[2]];
frachi = comm->zsplit[comm->myloc[2]+1];
partition_grid(nz,fraclo,frachi,shift,inzlo,inzhi);
} else {
fraclo = comm->mysplit[0][0];
frachi = comm->mysplit[0][1];
partition_grid(nx,fraclo,frachi,shift,inxlo,inxhi);
fraclo = comm->mysplit[1][0];
frachi = comm->mysplit[1][1];
partition_grid(ny,fraclo,frachi,shift,inylo,inyhi);
fraclo = comm->mysplit[2][0];
frachi = comm->mysplit[2][1];
partition_grid(nz,fraclo,frachi,shift,inzlo,inzhi);
}
ghost_grid(maxdist,extra);
// error check on size of grid stored by this proc
bigint total = (bigint)
(outxhi - outxlo + 1) * (outyhi - outylo + 1) * (outzhi - outzlo + 1);
if (total > MAXSMALLINT)
error->one(FLERR, "Too many owned+ghost grid3d points");
// default = caller grid is allocated to ghost grid
// used when computing pack/unpack lists in indices()
// NOTE: allow caller to override this
fullxlo = outxlo;
fullxhi = outxhi;
fullylo = outylo;
fullyhi = outyhi;
fullzlo = outzlo;
fullzhi = outzhi;
// initialize data structs
nswap = maxswap = 0;
swap = nullptr;
nsend = nrecv = ncopy = 0;
send = nullptr;
recv = nullptr;
copy = nullptr;
requests = nullptr;
requests_remap = nullptr;
xsplit = ysplit = zsplit = nullptr;
grid2proc = nullptr;
rcbinfo = nullptr;
nsend_remap = nrecv_remap = self_remap = 0;
send_remap = nullptr;
recv_remap = nullptr;
// 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();
// 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;
maxdist = 0.0;
stencil_grid_lo = stencil_grid_hi = 0;
stencil_atom_lo = stencil_atom_hi = 0;
shift_grid = shift_atom = 0.0;
zfactor = 1.0;
}
/* ----------------------------------------------------------------------
@ -367,7 +264,78 @@ void Grid3d::store(int ixlo, int ixhi, int iylo, int iyhi,
}
// ----------------------------------------------------------------------
// access Grid parameters
// set Grid parameters
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
maxdist = max distance outside proc subdomain a particle can be
------------------------------------------------------------------------- */
void Grid3d::set_distance(double distance)
{
maxdist = 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
------------------------------------------------------------------------- */
void Grid3d::set_stencil_grid(int lo, int hi)
{
stencil_grid_lo = lo;
stencil_grid_hi = hi;
}
/* ----------------------------------------------------------------------
# of grid cells beyond a particle's grid cell that caller updates
e.g. for smearing a point charge to the grid
can be different in lo vs hi direction
------------------------------------------------------------------------- */
void Grid3d::set_stencil_atom(int lo, int hi)
{
stencil_atom_lo = lo;
stencil_atom_hi = 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
------------------------------------------------------------------------- */
void Grid3d::set_shift_grid(double shift)
{
shift_grid = shift;
}
/* ----------------------------------------------------------------------
shift_atom = offset of 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
------------------------------------------------------------------------- */
void Grid3d::set_shift_atom(double shift)
{
shift_atom = shift;
}
/* ----------------------------------------------------------------------
zfactor > 1.0 when grid extends beyond simulation box bounds in Z
used by KSpace PPPM for 2d periodic slab geometries
default = 1.0 (no extension in Z)
------------------------------------------------------------------------- */
void Grid3d::set_zfactor(double factor)
{
zfactor = factor;
}
// ----------------------------------------------------------------------
// retrieve Grid parameters
// ----------------------------------------------------------------------
int Grid3d::identical(Grid3d *grid2)
@ -375,7 +343,7 @@ int Grid3d::identical(Grid3d *grid2)
int inxlo2,inxhi2,inylo2,inyhi2,inzlo2,inzhi2;
int outxlo2,outxhi2,outylo2,outyhi2,outzlo2,outzhi2;
grid2->get_bounds(inxlo2,inxhi2,inylo2,inyhi2,inzlo2,inzhi2);
grid2->get_bounds_owned(inxlo2,inxhi2,inylo2,inyhi2,inzlo2,inzhi2);
grid2->get_bounds_ghost(outxlo2,outxhi2,outylo2,outyhi2,outzlo2,outzhi2);
int flag = 0;
@ -404,8 +372,8 @@ void Grid3d::get_size(int &nxgrid, int &nygrid, int &nzgrid)
/* ---------------------------------------------------------------------- */
void Grid3d::get_bounds(int &xlo, int &xhi, int &ylo, int &yhi,
int &zlo, int &zhi)
void Grid3d::get_bounds_owned(int &xlo, int &xhi, int &ylo, int &yhi,
int &zlo, int &zhi)
{
xlo = inxlo;
xhi = inxhi;
@ -433,6 +401,113 @@ void Grid3d::get_bounds_ghost(int &xlo, int &xhi, int &ylo, int &yhi,
// also store comm and grid partitioning info
// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
setup grid partition for each proc = owned + ghost cells
return:
i xyz lohi = portion of global grid this proc owns, 0 <= index < N
o xyz lohi = owned + ghost grid cells in all directions
for periodic dims, o indices can be < 0 or >= N
for non-periodic dims, o indices will be >= 0 and < N
since no grid comm is done across non-periodic boundaries
------------------------------------------------------------------------- */
void Grid3d::setup_grid(int &ixlo, int &ixhi, int &iylo, int &iyhi,
int &izlo, int &izhi,
int &oxlo, int &oxhi, int &oylo, int &oyhi,
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);
fraclo = comm->ysplit[comm->myloc[1]];
frachi = comm->ysplit[comm->myloc[1]+1];
partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi);
fraclo = comm->zsplit[comm->myloc[2]];
frachi = comm->zsplit[comm->myloc[2]+1];
partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi);
} else {
fraclo = comm->mysplit[0][0];
frachi = comm->mysplit[0][1];
partition_grid(nx,fraclo,frachi,shift_grid,inxlo,inxhi);
fraclo = comm->mysplit[1][0];
frachi = comm->mysplit[1][1];
partition_grid(ny,fraclo,frachi,shift_grid,inylo,inyhi);
fraclo = comm->mysplit[2][0];
frachi = comm->mysplit[2][1];
partition_grid(nz,fraclo,frachi,shift_grid,inzlo,inzhi);
}
// extend owned grid bounds with ghost grid cells in each direction
ghost_grid();
// error check on size of grid stored by this proc
bigint total = (bigint)
(outxhi - outxlo + 1) * (outyhi - outylo + 1) * (outzhi - outzlo + 1);
if (total > MAXSMALLINT)
error->one(FLERR, "Too many owned+ghost grid3d points");
// default = caller grid is allocated to ghost grid
// used when computing pack/unpack lists in indices()
// NOTE: allow caller to override this
fullxlo = outxlo;
fullxhi = outxhi;
fullylo = outylo;
fullyhi = outyhi;
fullzlo = outzlo;
fullzhi = outzhi;
// initialize data structs
nswap = maxswap = 0;
swap = nullptr;
nsend = nrecv = ncopy = 0;
send = nullptr;
recv = nullptr;
copy = nullptr;
requests = nullptr;
requests_remap = nullptr;
xsplit = ysplit = zsplit = nullptr;
grid2proc = nullptr;
rcbinfo = nullptr;
nsend_remap = nrecv_remap = self_remap = 0;
send_remap = nullptr;
recv_remap = nullptr;
// 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();
// 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;
}
/* ----------------------------------------------------------------------
partition a global regular grid into one brick-shaped sub-grid per proc
if grid point is inside my sub-domain I own it,
@ -468,7 +543,7 @@ void Grid3d::partition_grid(int ngrid, double fraclo, double frachi,
ghost cell indices for periodic systems can be < 0 or >= N
------------------------------------------------------------------------- */
void Grid3d::ghost_grid(double maxdist, int extra)
void Grid3d::ghost_grid()
{
double *prd,*boxlo,*sublo,*subhi;
int triclinic = domain->triclinic;
@ -491,10 +566,11 @@ void Grid3d::ghost_grid(double maxdist, int extra)
if (triclinic == 0) dist[0] = dist[1] = dist[2] = maxdist;
else MathExtra::tribbox(domain->h,maxdist,&dist[0]);
// nlo,nhi = min/max index of global grid cell my owned atoms can be mapped to
// including up to maxdist displacement outside my subdomain
// extra = additional ghost layers required by called (e.g. finite diff stencil)
// max of the two quantities = ghost cell layers needed in each dim/dir
// lo/hi = min/max index of global grid cells my owned atoms can be mapped to
// includes effects of maxdist and shift_atom settings
// lo/hi can be further extended by stencil_atom and stencil_grid settings
// all those settings are set by caller
// ghost cell layers needed in each dim/dir = max of two extension effects
// OFFSET allows generation of negative indices with static_cast
// out xyz lo/hi = index range of owned + ghost cells
@ -503,22 +579,22 @@ void Grid3d::ghost_grid(double maxdist, int extra)
double dzinv = nz / prd[2];
int lo, hi;
lo = static_cast<int>((sublo[0]-dist[0]-boxlo[0]) * dxinv + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[0]+dist[0]-boxlo[0]) * dxinv + OFFSET) - OFFSET;
outxlo = MIN(lo, inxlo - extra);
outxhi = MAX(hi, inxhi + extra);
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;
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 + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[1]+dist[1]-boxlo[1]) * dyinv + OFFSET) - OFFSET;
outylo = MIN(lo, inylo - extra);
outyhi = MAX(hi, inyhi + extra);
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;
outxlo = MIN(lo - stencil_atom_lo, inylo - stencil_grid_lo);
outxhi = MAX(hi + stencil_atom_hi, inyhi + stencil_grid_hi);
// NOTE: need to account for zfactor here
lo = static_cast<int>((sublo[2]-dist[2]-boxlo[2]) * dzinv + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[2]+dist[2]-boxlo[2]) * dzinv + OFFSET) - OFFSET;
outzlo = MIN(lo, inzlo - extra);
outzhi = MAX(hi, inzhi + extra);
lo = static_cast<int>((sublo[2]-dist[2]-boxlo[2]) * dzinv + shift_atom + OFFSET) - OFFSET;
hi = static_cast<int>((subhi[2]+dist[2]-boxlo[2]) * dzinv + shift_atom + OFFSET) - OFFSET;
outxlo = MIN(lo - stencil_atom_lo, inzlo - stencil_grid_lo);
outxhi = MAX(hi + stencil_atom_hi, inzhi + stencil_grid_hi);
// limit out xyz lo/hi indices to global grid for non-periodic dims
@ -1427,7 +1503,8 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
// overlap_old = vector of overlap info in Overlap data struct
int oldbox[6];
old->get_bounds(oldbox[0],oldbox[1],oldbox[2],oldbox[3],oldbox[4],oldbox[5]);
old->get_bounds_owned(oldbox[0],oldbox[1],oldbox[2],oldbox[3],
oldbox[4],oldbox[5]);
pbc[0] = pbc[1] = pbc[2] = 0;
Overlap *overlap_old;
@ -1466,7 +1543,7 @@ void Grid3d::setup_remap(Grid3d *old, int &nremap_buf1, int &nremap_buf2)
// overlap_new = vector of overlap info in Overlap data struct
int newbox[6];
get_bounds(newbox[0],newbox[1],newbox[2],newbox[3],newbox[4],newbox[5]);
get_bounds_owned(newbox[0],newbox[1],newbox[2],newbox[3],newbox[4],newbox[5]);
pbc[0] = pbc[1] = pbc[2] = 0;
Overlap *overlap_new;
@ -1762,9 +1839,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,obox[0],obox[1]);
partition_grid(ny,ysplit[j],ysplit[j+1],shift,obox[2],obox[3]);
partition_grid(nz,zsplit[k],zsplit[k+1],shift,obox[4],obox[5]);
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]);
if (noverlap_list == maxoverlap_list) grow_overlap();
overlap_list[noverlap_list].proc = grid2proc[i][j][k];

37
src/grid3d.h
View File

@ -22,21 +22,31 @@ class Grid3d : protected Pointers {
public:
enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes
Grid3d(class LAMMPS *, MPI_Comm, int, int, int, double, int, double,
int &, int &, int &, int &, int &, int &,
int &, int &, int &, int &, int &, int &);
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(class LAMMPS *, MPI_Comm, int, int, int, int,
int, int, int, int, int, int, int, int, int, int, int, int,
int, int, int, int, int, int);
~Grid3d() 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_zfactor(double);
int identical(Grid3d *);
void get_size(int &, int &, int &);
void get_bounds(int &, int &, int &, int &, int &, int &);
void get_bounds_owned(int &, int &, int &, int &, int &, int &);
void get_bounds_ghost(int &, int &, int &, int &, int &, int &);
void setup_grid(int &, int &, int &, int &, int &, int &,
int &, int &, int &, int &, int &, int &);
void setup_comm(int &, int &);
int ghost_adjacent();
void forward_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
@ -54,9 +64,21 @@ class Grid3d : protected Pointers {
MPI_Comm gridcomm; // communicator for this class
// usually world, but MSM calls with subset
// inputs from caller via constructor
// inputs from caller
int nx, ny, nz; // size of global grid in all 3 dims
double maxdist; // distance owned atoms can move outside subdomain
int stencil_atom_lo,stencil_atom_hi; // grid cells accessed beyond atom's cell
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
// only affects extent of ghost cells
double zfactor; // multiplier on extend of grid in Z direction
// 1.0 = no extra grid cells in Z
// extent of my owned and ghost cells
int inxlo, inxhi; // inclusive extent of my grid chunk, 0 <= in <= N-1
int inylo, inyhi;
int inzlo, inzhi;
@ -67,9 +89,6 @@ class Grid3d : protected Pointers {
int fullylo, fullyhi; // can be same as out indices or larger
int fullzlo, fullzhi;
double shift; // location of grid point within grid cell
// only affects which proc owns grid cell
// -------------------------------------------
// internal variables for BRICK layout
// -------------------------------------------
@ -215,7 +234,7 @@ class Grid3d : protected Pointers {
// -------------------------------------------
void partition_grid(int, double, double, double, int &, int &);
void ghost_grid(double, int);
void ghost_grid();
void extract_comm_info();
void store(int, int, int, int, int, int, int, int, int, int, int, int,