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

more work on grid classes

Browse Source
This commit is contained in:
Steve Plimpton
2022-07-25 13:56:28 -06:00
parent 4a7daaa5d1
commit dff206f297
10 changed files with 924 additions and 237 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -39,7 +39,15 @@ using namespace FixConst;
static constexpr int MAXLINE = 256;
static constexpr int CHUNK = 1024;
// OFFSET avoids outside-of-box atoms being rounded to grid pts incorrectly
// SHIFT = 0.0 assigns atoms to lower-left grid pt
// SHIFT = 0.5 assigns atoms to nearest grid pt
// use SHIFT = 0.0 for now since it allows fix ave/chunk
// to spatially average consistent with the TTM grid
static constexpr int OFFSET = 16384;
static constexpr double SHIFT = 0.5;
/* ---------------------------------------------------------------------- */
@ -48,6 +56,11 @@ FixTTMGrid::FixTTMGrid(LAMMPS *lmp, int narg, char **arg) :
{
pergrid_flag = 1;
/*
if (outfile) error->all(FLERR,"Fix ttm/grid does not support outfile option - "
"use dump grid instead");
*/
skin_original = neighbor->skin;
}
@ -88,7 +101,8 @@ void FixTTMGrid::post_constructor()
if (infile) {
read_electron_temperatures(infile);
gc->forward_comm(Grid3d::FIX,this,1,sizeof(double),0,gc_buf1,gc_buf2,MPI_DOUBLE);
grid->forward_comm(Grid3d::FIX,this,1,sizeof(double),0,
grid_buf1,grid_buf2,MPI_DOUBLE);
}
}
@ -195,8 +209,8 @@ void FixTTMGrid::end_of_step()
flangevin[i][2]*v[i][2]);
}
gc->reverse_comm(Grid3d::FIX,this,1,sizeof(double),0,
gc_buf1,gc_buf2,MPI_DOUBLE);
grid->reverse_comm(Grid3d::FIX,this,1,sizeof(double),0,
grid_buf1,grid_buf2,MPI_DOUBLE);
// clang-format off
@ -248,7 +262,8 @@ void FixTTMGrid::end_of_step()
// communicate new T_electron values to ghost grid points
gc->forward_comm(Grid3d::FIX,this,1,sizeof(double),0,gc_buf1,gc_buf2,MPI_DOUBLE);
grid->forward_comm(Grid3d::FIX,this,1,sizeof(double),0,
grid_buf1,grid_buf2,MPI_DOUBLE);
}
// clang-format on
@ -365,7 +380,7 @@ void FixTTMGrid::write_electron_temperatures(const std::string &filename)
style);
}
gc->gather(Grid3d::FIX, this, 1, sizeof(double), 1, nullptr, MPI_DOUBLE);
grid->gather(Grid3d::FIX, this, 1, sizeof(double), 1, nullptr, MPI_DOUBLE);
if (comm->me == 0) fclose(FPout);
}
@ -424,42 +439,13 @@ void FixTTMGrid::unpack_reverse_grid(int /*flag*/, void *vbuf, int nlist, int *l
void FixTTMGrid::allocate_grid()
{
// 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
double maxdist = 0.5 * neighbor->skin;
comm->partition_grid(nxgrid, nygrid, nzgrid, 0.0, nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in,
nzhi_in);
grid = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, 1, maxdist, SHIFT,
nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in, nzhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
// nlo,nhi = min/max index of global grid pt my owned atoms can be mapped to
// finite difference stencil requires extra grid pt around my owned grid pts
// max of these 2 quantities is the ghost cells needed in each dim
// nlo_out,nhi_out = nlo_in,nhi_in + ghost cells
double *boxlo = domain->boxlo;
double *sublo = domain->sublo;
double *subhi = domain->subhi;
double dxinv = nxgrid / domain->xprd;
double dyinv = nxgrid / domain->yprd;
double dzinv = nxgrid / domain->zprd;
int nlo, nhi;
double cuthalf = 0.5 * neighbor->skin;
nlo = static_cast<int>((sublo[0] - cuthalf - boxlo[0]) * dxinv + shift) - OFFSET;
nhi = static_cast<int>((subhi[0] + cuthalf - boxlo[0]) * dxinv + shift) - OFFSET;
nxlo_out = MIN(nlo, nxlo_in - 1);
nxhi_out = MAX(nhi, nxhi_in + 1);
nlo = static_cast<int>((sublo[1] - cuthalf - boxlo[1]) * dyinv + shift) - OFFSET;
nhi = static_cast<int>((subhi[1] + cuthalf - boxlo[1]) * dyinv + shift) - OFFSET;
nylo_out = MIN(nlo, nylo_in - 1);
nyhi_out = MAX(nhi, nyhi_in + 1);
nlo = static_cast<int>((sublo[2] - cuthalf - boxlo[2]) * dzinv + shift) - OFFSET;
nhi = static_cast<int>((subhi[2] + cuthalf - boxlo[2]) * dzinv + shift) - OFFSET;
nzlo_out = MIN(nlo, nzlo_in - 1);
nzhi_out = MAX(nhi, nzhi_in + 1);
// set ngridout and ngridmine and error check
bigint totalmine;
totalmine =
@ -470,13 +456,12 @@ void FixTTMGrid::allocate_grid()
totalmine = (bigint) (nxhi_in - nxlo_in + 1) * (nyhi_in - nylo_in + 1) * (nzhi_in - nzlo_in + 1);
ngridmine = totalmine;
gc = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, nxlo_in, nxhi_in, nylo_in, nyhi_in, nzlo_in,
nzhi_in, nxlo_out, nxhi_out, nylo_out, nyhi_out, nzlo_out, nzhi_out);
// setup grid communication and allocate grid data structs
gc->setup(ngc_buf1, ngc_buf2);
grid->setup(ngrid_buf1, ngrid_buf2);
memory->create(gc_buf1, ngc_buf1, "ttm/grid:gc_buf1");
memory->create(gc_buf2, ngc_buf2, "ttm/grid:gc_buf2");
memory->create(grid_buf1, ngrid_buf1, "ttm/grid:grid_buf1");
memory->create(grid_buf2, ngrid_buf2, "ttm/grid:grid_buf2");
memory->create3d_offset(T_electron_old, nzlo_out, nzhi_out, nylo_out, nyhi_out, nxlo_out,
nxhi_out, "ttm/grid:T_electron_old");
@ -492,9 +477,9 @@ void FixTTMGrid::allocate_grid()
void FixTTMGrid::deallocate_grid()
{
delete gc;
memory->destroy(gc_buf1);
memory->destroy(gc_buf2);
delete grid;
memory->destroy(grid_buf1);
memory->destroy(grid_buf2);
memory->destroy3d_offset(T_electron_old, nzlo_out, nylo_out, nxlo_out);
memory->destroy3d_offset(T_electron, nzlo_out, nylo_out, nxlo_out);
@ -519,7 +504,7 @@ void FixTTMGrid::write_restart(FILE *fp)
// gather rest of rlist on proc 0 as global grid values
gc->gather(Grid3d::FIX, this, 1, sizeof(double), 0, &rlist[4], MPI_DOUBLE);
//grid->gather(Grid3d::FIX, this, 1, sizeof(double), 0, &rlist[4], MPI_DOUBLE);
if (comm->me == 0) {
int size = rsize * sizeof(double);
@ -570,7 +555,8 @@ void FixTTMGrid::restart(char *buf)
// communicate new T_electron values to ghost grid points
gc->forward_comm(Grid3d::FIX, this, 1, sizeof(double), 0, gc_buf1, gc_buf2, MPI_DOUBLE);
grid->forward_comm(Grid3d::FIX, this, 1, sizeof(double), 0,
grid_buf1, grid_buf2, MPI_DOUBLE);
}
/* ----------------------------------------------------------------------
@ -650,7 +636,7 @@ int FixTTMGrid::get_grid_by_name(char *name, int &dim)
void *FixTTMGrid::get_grid_by_index(int index)
{
if (index == 0) return gc;
if (index == 0) return grid;
return nullptr;
}

6
src/EXTRA-FIX/fix_ttm_grid.h
View File

@ -61,9 +61,9 @@ class FixTTMGrid : public FixTTM {
double skin_original;
FILE *FPout;
class Grid3d *gc;
int ngc_buf1, ngc_buf2;
double *gc_buf1, *gc_buf2;
class Grid3d *grid;
int ngrid_buf1, ngrid_buf2;
double *grid_buf1, *grid_buf2;
void allocate_grid() override;
void deallocate_grid() override;

37
src/comm.cpp
View File

@ -806,6 +806,7 @@ int Comm::coord2proc(double *x, int &igx, int &igy, int &igz)
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
use nz=1 for 2d grid
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
@ -859,42 +860,6 @@ void Comm::partition_grid(int nx, int ny, int nz, double zfactor,
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;
}
}
*/
}
/* ----------------------------------------------------------------------

423
src/compute_property_grid.cpp
View File

@ -31,7 +31,9 @@ enum { ID, X, Y, Z, XS, YS, ZS, XC, YC, ZC, XSC, YSC, ZSC };
/* ---------------------------------------------------------------------- */
ComputePropertyGrid::ComputePropertyGrid(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), pack_choice(nullptr)
Compute(lmp, narg, arg), pack_choice(nullptr),
grid2d(nullptr), grid3d(nullptr),
vec2d(nullptr), array2d(nullptr), vec3d(nullptr), array3d(nullptr)
{
if (narg < 7) error->all(FLERR, "Illegal compute property/grid command");
@ -53,53 +55,78 @@ ComputePropertyGrid::ComputePropertyGrid(LAMMPS *lmp, int narg, char **arg) :
pack_choice = new FnPtrPack[nvalues];
for (int iarg = 6; iarg < narg; iarg++) {
int jarg = iarg - 6;
if (strcmp(arg[iarg], "id") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_id;
pack_choice[jarg] = &ComputePropertyGrid::pack_id;
} else if (strcmp(arg[iarg], "ix") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_ix;
} else if (strcmp(arg[iarg], "iy") == 0) {
pack_choice[jarg] = &ComputePropertyGrid::pack_iy;
} else if (strcmp(arg[iarg], "iz") == 0) {
if (dimension == 2)
error->all(FLERR,"Compute property/grid for 2d cannot use z coord");
pack_choice[jarg] = &ComputePropertyGrid::pack_iz;
} else if (strcmp(arg[iarg], "x") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_x;
pack_choice[jarg] = &ComputePropertyGrid::pack_x;
} else if (strcmp(arg[iarg], "y") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_y;
pack_choice[jarg] = &ComputePropertyGrid::pack_y;
} else if (strcmp(arg[iarg], "z") == 0) {
if (dimension == 2)
error->all(FLERR,"Compute property/grid for 2d cannot use z coord");
pack_choice[iarg] = &ComputePropertyGrid::pack_z;
pack_choice[jarg] = &ComputePropertyGrid::pack_z;
} else if (strcmp(arg[iarg], "xs") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_xs;
pack_choice[jarg] = &ComputePropertyGrid::pack_xs;
} else if (strcmp(arg[iarg], "ys") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_ys;
pack_choice[jarg] = &ComputePropertyGrid::pack_ys;
} else if (strcmp(arg[iarg], "zs") == 0) {
if (dimension == 2)
error->all(FLERR,"Compute property/grid for 2d cannot use z coord");
pack_choice[iarg] = &ComputePropertyGrid::pack_zs;
pack_choice[jarg] = &ComputePropertyGrid::pack_zs;
} else if (strcmp(arg[iarg], "xc") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_xc;
pack_choice[jarg] = &ComputePropertyGrid::pack_xc;
} else if (strcmp(arg[iarg], "yc") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_yc;
pack_choice[jarg] = &ComputePropertyGrid::pack_yc;
} else if (strcmp(arg[iarg], "zc") == 0) {
if (dimension == 2)
error->all(FLERR,"Compute property/grid for 2d cannot use z coord");
pack_choice[iarg] = &ComputePropertyGrid::pack_zc;
pack_choice[jarg] = &ComputePropertyGrid::pack_zc;
} else if (strcmp(arg[iarg], "xsc") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_xsc;
pack_choice[jarg] = &ComputePropertyGrid::pack_xsc;
} else if (strcmp(arg[iarg], "ysc") == 0) {
pack_choice[iarg] = &ComputePropertyGrid::pack_ysc;
pack_choice[jarg] = &ComputePropertyGrid::pack_ysc;
} else if (strcmp(arg[iarg], "zsc") == 0) {
if (dimension == 2)
error->all(FLERR,"Compute property/grid for 2d cannot use z coord");
pack_choice[iarg] = &ComputePropertyGrid::pack_zsc;
pack_choice[jarg] = &ComputePropertyGrid::pack_zsc;
} else error->all(FLERR, "Illegal compute property/grid command");
}
// instantiate the Grid class and allocate per-grid memory
// NOTE: need new memory create methods for 2d
if (dimension == 2) {
// instanttiate the Grid2d/3d class
} else {
grid3d = new Grid3d(lmp, world, nx, ny, nz, 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);
if (nvalues == 1)
memory->create3d_offset(vec3d, nzlo_out, nzhi_out, nylo_out,
nyhi_out, nxlo_out,
nxhi_out, "property/grid:vec3d");
else
memory->create4d_offset_last(array3d, nzlo_out, nzhi_out, nylo_out,
nyhi_out, nxlo_out,
nxhi_out, nvalues, "property/grid:array3d");
}
}
/* ---------------------------------------------------------------------- */
@ -107,6 +134,13 @@ ComputePropertyGrid::ComputePropertyGrid(LAMMPS *lmp, int narg, char **arg) :
ComputePropertyGrid::~ComputePropertyGrid()
{
delete[] pack_choice;
delete grid2d;
delete grid3d;
//memory->destroy2d_offset(vec2d);
//memory->destroy2d_offset(array2d);
memory->destroy3d_offset(vec3d,nzlo_out,nylo_out,nxlo_out);
memory->destroy4d_offset_last(array3d,nzlo_out,nylo_out,nxlo_out);
}
/* ---------------------------------------------------------------------- */
@ -119,9 +153,9 @@ void ComputePropertyGrid::compute_pergrid()
// may change between compute invocations due to load balancing
if (dimension == 2)
grid2d->query_in_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in);
grid2d->query_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in);
else
grid3d->query_in_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
grid3d->query_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
// reallocate data vector or array if changed
@ -179,7 +213,8 @@ void *ComputePropertyGrid::get_grid_by_index(int index)
int ComputePropertyGrid::get_griddata_by_name(int igrid, char *name, int &ncol)
{
if (igrid == 0 && strcmp(name,"data") == 0) {
ncol = nvalues;
if (nvalues == 1) ncol = 0;
else ncol = nvalues;
return 0;
}
@ -195,10 +230,10 @@ void *ComputePropertyGrid::get_griddata_by_index(int index)
{
if (index == 0) {
if (dimension == 2) {
if (nvalues == 0) return vec2d;
if (nvalues == 1) return vec2d;
else return array2d;
} else {
if (nvalues == 0) return vec3d;
if (nvalues == 1) return vec3d;
else return array3d;
}
}
@ -219,30 +254,27 @@ double ComputePropertyGrid::memory_usage()
/* ----------------------------------------------------------------------
one method for every keyword compute property/grid can output
packed into buf starting at n with stride nvalues
------------------------------------------------------------------------- */
void ComputePropertyGrid::pack_id(int n)
{
if (dimension == 2) {
if (n == 0) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = iy*nx + ix + 1;
} else {
n--;
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = iy*nx + ix + 1;
}
} else if (dimension == 3) {
if (n == 0) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = iz*ny*nx + iy*nx + ix + 1;
} else {
n--;
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
@ -253,72 +285,409 @@ void ComputePropertyGrid::pack_id(int n)
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_ix(int n)
{
if (dimension == 2) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = ix + 1;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = ix + 1;
}
} else if (dimension == 3) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = ix + 1;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = ix + 1;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_iy(int n)
{
if (dimension == 2) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = iy + 1;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = iy + 1;
}
} else if (dimension == 3) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = iy + 1;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = iy + 1;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_iz(int n)
{
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = iz + 1;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = iz + 1;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_x(int n)
{
double boxlo,dx;
if (dimension == 2) {
grid2d->query_box(0,boxlo,dx);
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = boxlo + ix*dx;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = boxlo + ix*dx;
}
} else if (dimension == 3) {
grid3d->query_box(0,boxlo,dx);
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = boxlo + ix*dx;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = boxlo + ix*dx;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_y(int n)
{
double boxlo,dy;
if (dimension == 2) {
grid2d->query_box(1,boxlo,dy);
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = boxlo + iy*dy;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = boxlo + iy*dy;
}
} else if (dimension == 3) {
grid3d->query_box(1,boxlo,dy);
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = boxlo + iy*dy;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = boxlo + iy*dy;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_z(int n)
{
double boxlo,dz;
grid3d->query_box(2,boxlo,dz);
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = boxlo + iz*dz;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = boxlo + iz*dz;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_xs(int n)
{
double dx = 1.0/nx;
if (dimension == 2) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = ix*dx;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = ix*dx;
}
} else if (dimension == 3) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = ix*dx;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = ix*dx;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_ys(int n)
{
double dy = 1.0/ny;
if (dimension == 2) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = iy*dy;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = iy*dy;
}
} else if (dimension == 3) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = iy*dy;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = iy*dy;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_zs(int n)
{
double dz = 1.0/nz;
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = iz*dz;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = iz*dz;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_xc(int n)
{
double boxlo,dx;
if (dimension == 2) {
grid2d->query_box(0,boxlo,dx);
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = boxlo + (ix+0.5)*dx;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = boxlo + (ix+0.5)*dx;
}
} else if (dimension == 3) {
grid3d->query_box(0,boxlo,dx);
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = boxlo + (ix+0.5)*dx;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = boxlo + (ix+0.5)*dx;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_yc(int n)
{
double boxlo,dy;
if (dimension == 2) {
grid2d->query_box(1,boxlo,dy);
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = boxlo + (iy+0.5)*dy;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = boxlo + (iy+0.5)*dy;
}
} else if (dimension == 3) {
grid3d->query_box(1,boxlo,dy);
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = boxlo + (iy+0.5)*dy;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = boxlo + (iy+0.5)*dy;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_zc(int n)
{
double boxlo,dz;
grid3d->query_box(2,boxlo,dz);
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = boxlo + (iz+0.5)*dz;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = boxlo + (iz+0.5)*dz;
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_xsc(int n)
{
double dx = 1.0/nx;
if (dimension == 2) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = (ix+0.5)*dx;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = (ix+0.5)*dx;
}
} else if (dimension == 3) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = (ix+0.5)*dx;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = (ix+0.5)*dx;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_ysc(int n)
{
double dy = 1.0/ny;
if (dimension == 2) {
if (nvalues == 0) {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] = (iy+0.5)*dy;
} else {
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array2d[iy][ix][n] = (iy+0.5)*dy;
}
} else if (dimension == 3) {
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = (iy+0.5)*dy;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = (iy+0.5)*dy;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputePropertyGrid::pack_zsc(int n)
{
double dz = 1.0/nz;
if (nvalues == 0) {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] = (iz+0.5)*dz;
} else {
for (int iz = nzlo_in; iz <= nzhi_in; iz++)
for (int iy = nylo_in; iy <= nyhi_in; iy++)
for (int ix = nxlo_in; ix <= nxhi_in; ix++)
array3d[iz][iy][ix][n] = (iz+0.5)*dz;
}
}

5
src/compute_property_grid.h
View File

@ -47,6 +47,7 @@ class ComputePropertyGrid : public Compute {
class Grid3d *grid3d;
int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
double **vec2d,***vec3d;
double ***array2d,****array3d;
@ -56,6 +57,10 @@ class ComputePropertyGrid : public Compute {
void pack_id(int);
void pack_ix(int);
void pack_iy(int);
void pack_iz(int);
void pack_x(int);
void pack_y(int);
void pack_z(int);

14
src/dump_grid.cpp
View File

@ -274,7 +274,7 @@ void DumpGrid::init_style()
int nx,ny,nz,nxtmp,nytmp,nztmp;
for (int i = 1; i < nfield; i++) {
for (int i = 0; i < nfield; i++) {
if (dimension == 2) {
if (field2source[i] == COMPUTE) {
icompute = compute[field2index[i]];
@ -283,9 +283,9 @@ void DumpGrid::init_style()
ifix = fix[field2index[i]];
grid2d = (Grid2d *) ifix->get_grid_by_index(field2grid[i]);
}
if (i == 0) grid2d->query_global_size(nx,ny);
if (i == 0) grid2d->query_size(nx,ny);
else {
grid2d->query_global_size(nxtmp,nytmp);
grid2d->query_size(nxtmp,nytmp);
if (nxtmp != nx || nytmp != ny)
error->all(FLERR,"Dump grid field grid sizes do not match");
}
@ -298,9 +298,9 @@ void DumpGrid::init_style()
ifix = fix[field2index[i]];
grid3d = (Grid3d *) ifix->get_grid_by_index(field2grid[i]);
}
if (i == 0) grid3d->query_global_size(nx,ny,nz);
if (i == 0) grid3d->query_size(nx,ny,nz);
else {
grid3d->query_global_size(nxtmp,nytmp,nztmp);
grid3d->query_size(nxtmp,nytmp,nztmp);
if (nxtmp != nx || nytmp != ny || nztmp != nz)
error->all(FLERR,"Dump grid field grid sizes do not match");
}
@ -535,7 +535,7 @@ int DumpGrid::count()
else if (field2source[0] == FIX)
grid2d = (Grid2d *)
fix[field2index[0]]->get_grid_by_index(field2grid[0]);
grid2d->query_in_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in);
grid2d->query_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in);
} else {
if (field2source[0] == COMPUTE)
grid3d = (Grid3d *)
@ -543,7 +543,7 @@ int DumpGrid::count()
else if (field2source[0] == FIX)
grid3d = (Grid3d *)
fix[field2index[0]]->get_grid_by_index(field2grid[0]);
grid3d->query_in_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
grid3d->query_bounds(nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in);
}
// invoke Computes for per-grid quantities

288
src/grid2d.cpp
View File

@ -15,6 +15,7 @@
#include "grid2d.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "irregular.h"
#include "pair.h"
@ -28,6 +29,8 @@ enum{REGULAR,TILED};
#define DELTA 16
static constexpr int OFFSET = 16384;
/* ----------------------------------------------------------------------
NOTES
tiled implementation only currently works for RCB, not general tiled
@ -38,9 +41,108 @@ enum{REGULAR,TILED};
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
constructor called by all classes except MSM
constructor called by all classes except PPPM and MSM
gcomm = world communicator
gn xy = size of global grid
gnx, gny = 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 = 0.0 for grid pt in lower-left corner of grid cell, 0.5 for center
return:
i xy lohi = portion of global grid this proc owns, 0 <= index < N
o xy lohi = owned + ghost grid cells needed in all directions
for non-periodic dims, o indices will not be < 0 or >= N,
since no grid communication is done across non-periodic boundaries
------------------------------------------------------------------------- */
Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm,
int gnx, int gny,
double maxdist, int extra, double shift,
int &ixlo, int &ixhi, int &iylo, int &iyhi,
int &oxlo, int &oxhi, int &oylo, int &oyhi)
: Pointers(lmp)
{
// store commnicator and global grid size
// set layout mode
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
ngrid[0] = nx; ngrid[1] = ny;
if (comm->layout == Comm::LAYOUT_TILED) layout = TILED;
else layout = REGULAR;
// partition global grid across procs
// i xyz lo/hi = lower/upper bounds of global grid this proc owns
// indices range from 0 to N-1 inclusive in each dim
int tmp1,tmp2;
comm->partition_grid(nx, ny, 1, 0.0, ixlo, ixhi, iylo, iyhi, tmp1, tmp2);
// nlo,nhi = min/max index of global grid pt my owned atoms can be mapped to
// finite difference stencil requires extra grid pt around my owned grid pts
// max of these 2 quantities is the ghost cells needed in each dim
// o xyz lo/hi = owned + ghost cells
memcpy(boxlo,domain->boxlo,2*sizeof(double));
memcpy(prd,domain->prd,2*sizeof(double));
double *sublo = domain->sublo;
double *subhi = domain->subhi;
int *periodicity = domain->periodicity;
double dxinv = nx / domain->prd[0];
double dyinv = ny / domain->prd[1];
double SHIFT = OFFSET + shift;
int nlo, nhi;
nlo = static_cast<int>((sublo[0]-maxdist-boxlo[0]) * dxinv + SHIFT) - OFFSET;
nhi = static_cast<int>((subhi[0]+maxdist-boxlo[0]) * dxinv + SHIFT) - OFFSET;
oxlo = MIN(nlo, ixlo - extra);
oxhi = MAX(nhi, ixhi + extra);
nlo = static_cast<int>((sublo[1]-maxdist-boxlo[1]) * dyinv + SHIFT) - OFFSET;
nhi = static_cast<int>((subhi[1]+maxdist-boxlo[1]) * dyinv + SHIFT) - OFFSET;
oylo = MIN(nlo, iylo - extra);
oyhi = MAX(nhi, iyhi + extra);
// limit o xyz lo/hi indices for non-periodic dimensions
if (!periodicity[0]) {
oxlo = MAX(1,oxlo);
oxhi = MIN(gnx-1,oxhi);
}
if (!periodicity[1]) {
oylo = MAX(1,oylo);
oyhi = MIN(gnx-1,oyhi);
}
// store grid bounds and proc neighs
if (layout == REGULAR) {
int (*procneigh)[2] = comm->procneigh;
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1]);
} else {
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
0,0,0,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,
@ -54,21 +156,38 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm,
int oxlo, int oxhi, int oylo, int oyhi)
: Pointers(lmp)
{
// store commnicator and global grid size
// set layout mode
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
ngrid[0] = nx; ngrid[1] = ny;
if (comm->layout == Comm::LAYOUT_TILED) layout = TILED;
else layout = REGULAR;
memcpy(boxlo,domain->boxlo,2*sizeof(double));
memcpy(prd,domain->prd,2*sizeof(double));
// store grid bounds and proc neighs
if (layout == REGULAR) {
int (*procneigh)[2] = comm->procneigh;
initialize(gcomm,gnx,gny,
ixlo,ixhi,iylo,iyhi,oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1]);
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1]);
} else {
initialize(gcomm,gnx,gny,
ixlo,ixhi,iylo,iyhi,oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
0,0,0,0);
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
0,0,0,0);
}
}
@ -77,7 +196,7 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm,
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
gn xy = size of global grid
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
@ -91,34 +210,48 @@ Grid2d::Grid2d(LAMMPS *lmp, MPI_Comm gcomm, int flag,
int exlo, int exhi, int eylo, int eyhi)
: Pointers(lmp)
{
// store commnicator and global grid size
// set layout mode
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
ngrid[0] = nx; ngrid[1] = ny;
if (comm->layout == Comm::LAYOUT_TILED) layout = TILED;
else layout = REGULAR;
memcpy(boxlo,domain->boxlo,2*sizeof(double));
memcpy(prd,domain->prd,2*sizeof(double));
// store grid bounds and proc neighs
if (flag == 1) {
if (layout == REGULAR) {
// this assumes gcomm = world
int (*procneigh)[2] = comm->procneigh;
initialize(gcomm,gnx,gny,
ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1]);
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1]);
} else {
initialize(gcomm,gnx,gny,
ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi,
0,0,0,0);
store(ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi,
0,0,0,0);
}
} else if (flag == 2) {
if (layout == REGULAR) {
initialize(gcomm,gnx,gny,
ixlo,ixhi,iylo,iyhi,
oxlo,oxhi,oylo,oyhi,
oxlo,oxhi,oylo,oyhi,
exlo,exhi,eylo,eyhi);
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 with neighbor procs");
}
@ -160,20 +293,11 @@ Grid2d::~Grid2d()
store constructor args in local variables
------------------------------------------------------------------------- */
void Grid2d::initialize(MPI_Comm gcomm,
int gnx, int gny,
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)
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)
{
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
inxlo = ixlo;
inxhi = ixhi;
inylo = iylo;
@ -212,7 +336,7 @@ void Grid2d::initialize(MPI_Comm gcomm,
/* ---------------------------------------------------------------------- */
void Grid2d::query_global_size(int &nxgrid, int &nygrid)
void Grid2d::query_size(int &nxgrid, int &nygrid)
{
nxgrid = nx;
nygrid = ny;
@ -220,7 +344,7 @@ void Grid2d::query_global_size(int &nxgrid, int &nygrid)
/* ---------------------------------------------------------------------- */
void Grid2d::query_in_bounds(int &xlo, int &xhi, int &ylo, int &yhi)
void Grid2d::query_bounds(int &xlo, int &xhi, int &ylo, int &yhi)
{
xlo = inxlo;
xhi = inxhi;
@ -230,6 +354,14 @@ void Grid2d::query_in_bounds(int &xlo, int &xhi, int &ylo, int &yhi)
/* ---------------------------------------------------------------------- */
void Grid2d::query_box(int dim, double &lo, double &delta)
{
lo = boxlo[dim];
delta = prd[dim] / ngrid[dim];
}
/* ---------------------------------------------------------------------- */
void Grid2d::setup(int &nbuf1, int &nbuf2)
{
if (layout == REGULAR) setup_regular(nbuf1,nbuf2);
@ -1022,6 +1154,76 @@ reverse_comm_tiled(T *ptr, int nper, int nbyte, int which,
}
}
/* ----------------------------------------------------------------------
gather global grid values to proc 0, one grid chunk at a time
proc 0 pings each proc for its grid chunk
pack/unpack operations are performed by caller via callbacks
caller can decide whether to store chunks, output them, etc
------------------------------------------------------------------------- */
void Grid2d::gather(int /*caller*/, void *ptr, int nper, int nbyte,
int which, void *buf, MPI_Datatype datatype)
{
int me = comm->me;
Fix *fptr = (Fix *) ptr;
// maxsize = max grid data owned by any proc
int mysize = (inxhi-inxlo+1) * (inyhi-inylo+1);
mysize *= nper;
int maxsize;
MPI_Allreduce(&mysize,&maxsize,1,MPI_INT,MPI_MAX,world);
// pack my data via callback to caller
char *mybuf;
if (me == 0) memory->create(mybuf,maxsize*nbyte,"grid2d:mybuf");
else memory->create(mybuf,mysize*nbyte,"grid2d:mybuf");
fptr->pack_gather_grid(which,mybuf);
// ping each proc for its data
// unpack into full buffer via callback to caller
int xlo,xhi,ylo,yhi,zlo,zhi,tmp;
int bounds[4];
if (me == 0) {
MPI_Status status;
MPI_Request request;
for (int iproc = 0; iproc < nprocs; iproc++) {
if (iproc) {
MPI_Irecv(mybuf,maxsize,datatype,iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
MPI_Wait(&request,&status);
MPI_Recv(bounds,4,MPI_INT,iproc,0,world,&status);
xlo = bounds[0];
xhi = bounds[1];
ylo = bounds[2];
yhi = bounds[3];
} else {
xlo = inxlo;
xhi = inxhi;
ylo = inylo;
yhi = inyhi;
}
fptr->unpack_gather_grid(which,mybuf,buf,xlo,xhi,ylo,yhi,0,0);
}
} else {
MPI_Recv(&tmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
MPI_Rsend(mybuf,mysize,datatype,0,0,world);
bounds[0] = inxlo;
bounds[1] = inxhi;
bounds[2] = inylo;
bounds[3] = inyhi;
MPI_Send(bounds,4,MPI_INT,0,0,world);
}
memory->destroy(mybuf);
}
/* ----------------------------------------------------------------------
create swap stencil for grid own/ghost communication
swaps covers all 2 dimensions and both directions

25
src/grid2d.h
View File

@ -22,16 +22,21 @@ class Grid2d : protected Pointers {
public:
enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes
Grid2d(class LAMMPS *, MPI_Comm, int, int, double, int, double,
int &, int &, int &, int &,
int &, int &, 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 query_global_size(int &, int &);
void query_in_bounds(int &, int &, int &, int &);
void query_size(int &, int &);
void query_bounds(int &, int &, int &, int &);
void query_box(int, double &, double &);
void setup(int &, int &);
int ghost_adjacent();
void forward_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
void reverse_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
void gather(int, void *, int, int, int, void *, MPI_Datatype);
protected:
int me, nprocs;
@ -39,18 +44,20 @@ class Grid2d : protected Pointers {
MPI_Comm gridcomm; // communicator for this class
// usually world, but MSM calls with subset
int ngrid[2]; // global grid size
double boxlo[2]; // current box that grid is mapped to
double prd[2];
// inputs from caller via constructor
int nx, ny, nz; // size of global grid in all 3 dims
int nx, ny; // size of global grid in both dims
int inxlo, inxhi; // inclusive extent of my grid chunk
int inylo, inyhi; // 0 <= in <= N-1
int inzlo, inzhi;
int outxlo, outxhi; // inclusive extent of my grid chunk plus
int outylo, outyhi; // ghost cells in all 6 directions
int outzlo, outzhi; // lo indices can be < 0, hi indices can be >= N
int outylo, outyhi; // ghost cells in all 4 directions
// lo indices can be < 0, hi indices can be >= N
int fullxlo, fullxhi; // extent of grid chunk that caller stores
int fullylo, fullyhi; // can be same as out indices or larger
int fullzlo, fullzhi;
// -------------------------------------------
// internal variables for REGULAR layout
@ -172,8 +179,8 @@ class Grid2d : protected Pointers {
// internal methods
// -------------------------------------------
void initialize(MPI_Comm, int, int, int, int, int, int, int, int, int, int,
int, int, int, int, int, int, int, int);
void store(int, int, int, int, int, int, int, int,
int, int, int, int, int, int, int, int);
virtual void setup_regular(int &, int &);
virtual void setup_tiled(int &, int &);
void ghost_box_drop(int *, int *);

250
src/grid3d.cpp
View File

@ -15,6 +15,7 @@
#include "grid3d.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "irregular.h"
#include "pair.h"
@ -28,6 +29,8 @@ enum{REGULAR,TILED};
#define DELTA 16
static constexpr int OFFSET = 16384;
/* ----------------------------------------------------------------------
NOTES
tiled implementation only currently works for RCB, not general tiled
@ -38,9 +41,120 @@ enum{REGULAR,TILED};
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
constructor called by all classes except MSM
constructor called by all classes except PPPM and MSM
gcomm = world communicator
gn xyz = size of global grid
gnx, gny, gnz = 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 = 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 non-periodic dims, o indices will not be < 0 or >= N,
since no grid communication is done across non-periodic boundaries
------------------------------------------------------------------------- */
Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm,
int gnx, int gny, int gnz,
double maxdist, int extra, double shift,
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)
{
// store commnicator and global grid size
// set layout mode
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
nz = gnz;
ngrid[0] = nx; ngrid[1] = ny; ngrid[2] = nz;
if (comm->layout == Comm::LAYOUT_TILED) layout = TILED;
else layout = REGULAR;
// partition global grid across procs
// i xyz lo/hi = lower/upper bounds of global grid this proc owns
// indices range from 0 to N-1 inclusive in each dim
comm->partition_grid(nx, ny, nz, 0.0, ixlo, ixhi, iylo, iyhi, izlo, izhi);
// nlo,nhi = min/max index of global grid pt my owned atoms can be mapped to
// finite difference stencil requires extra grid pt around my owned grid pts
// max of these 2 quantities is the ghost cells needed in each dim
// o xyz lo/hi = owned + ghost cells
memcpy(boxlo,domain->boxlo,3*sizeof(double));
memcpy(prd,domain->prd,3*sizeof(double));
double *sublo = domain->sublo;
double *subhi = domain->subhi;
int *periodicity = domain->periodicity;
double dxinv = nx / prd[0];
double dyinv = ny / prd[1];
double dzinv = nz / prd[2];;
double SHIFT = OFFSET + shift;
int nlo, nhi;
nlo = static_cast<int>((sublo[0]-maxdist-boxlo[0]) * dxinv + SHIFT) - OFFSET;
nhi = static_cast<int>((subhi[0]+maxdist-boxlo[0]) * dxinv + SHIFT) - OFFSET;
oxlo = MIN(nlo, ixlo - extra);
oxhi = MAX(nhi, ixhi + extra);
nlo = static_cast<int>((sublo[1]-maxdist-boxlo[1]) * dyinv + SHIFT) - OFFSET;
nhi = static_cast<int>((subhi[1]+maxdist-boxlo[1]) * dyinv + SHIFT) - OFFSET;
oylo = MIN(nlo, iylo - extra);
oyhi = MAX(nhi, iyhi + extra);
nlo = static_cast<int>((sublo[2]-maxdist-boxlo[2]) * dzinv + SHIFT) - OFFSET;
nhi = static_cast<int>((subhi[2]+maxdist-boxlo[2]) * dzinv + SHIFT) - OFFSET;
ozlo = MIN(nlo, izlo - extra);
ozhi = MAX(nhi, izhi + extra);
// limit o xyz lo/hi indices for non-periodic dimensions
if (!periodicity[0]) {
oxlo = MAX(1,oxlo);
oxhi = MIN(gnx-1,oxhi);
}
if (!periodicity[1]) {
oylo = MAX(1,oylo);
oyhi = MIN(gnx-1,oyhi);
}
if (!periodicity[2]) {
ozlo = MAX(1,ozlo);
ozhi = MIN(gnx-1,ozhi);
}
// store grid bounds and proc neighs
if (layout == REGULAR) {
int (*procneigh)[2] = comm->procneigh;
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1],
procneigh[2][0],procneigh[2][1]);
} else {
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
0,0,0,0,0,0);
}
}
/* ----------------------------------------------------------------------
constructor called by PPPM classes
gcomm = world communicator
gnx, gny, gnz = size of global grid
i xyz lohi = portion of global grid this proc owns, 0 <= index < N
o xyz lohi = owned grid portion + ghost grid cells needed in all directions
if o indices are < 0 or hi indices are >= N,
@ -54,24 +168,40 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm,
int oxlo, int oxhi, int oylo, int oyhi, int ozlo, int ozhi)
: Pointers(lmp)
{
// store commnicator and global grid size
// set layout mode
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
nz = gnz;
ngrid[0] = nx; ngrid[1] = ny; ngrid[2] = nz;
if (comm->layout == Comm::LAYOUT_TILED) layout = TILED;
else layout = REGULAR;
memcpy(boxlo,domain->boxlo,3*sizeof(double));
memcpy(prd,domain->prd,3*sizeof(double));
// store grid bounds and proc neighs
if (layout == REGULAR) {
int (*procneigh)[2] = comm->procneigh;
initialize(gcomm,gnx,gny,gnz,
ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1],
procneigh[2][0],procneigh[2][1]);
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1],
procneigh[2][0],procneigh[2][1]);
} else {
initialize(gcomm,gnx,gny,gnz,
ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
0,0,0,0,0,0);
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
0,0,0,0,0,0);
}
}
@ -80,7 +210,7 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm,
gcomm = world communicator or sub-communicator for a hierarchical grid
flag = 1 if e xyz lohi values = larger grid stored by caller in gcomm = world
flag = 2 if e xyz lohi values = 6 neighbor procs in gcomm
gn xyz = size of global grid
gnx, gny, gnz = size of global grid
i xyz lohi = portion of global grid this proc owns, 0 <= index < N
o xyz lohi = owned grid portion + ghost grid cells needed in all directions
e xyz lohi for flag = 1: extent of larger grid stored by caller
@ -94,35 +224,50 @@ Grid3d::Grid3d(LAMMPS *lmp, MPI_Comm gcomm, int flag,
int exlo, int exhi, int eylo, int eyhi, int ezlo, int ezhi)
: Pointers(lmp)
{
// store commnicator and global grid size
// set layout mode
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
nz = gnz;
ngrid[0] = nx; ngrid[1] = ny; ngrid[2] = nz;
if (comm->layout == Comm::LAYOUT_TILED) layout = TILED;
else layout = REGULAR;
memcpy(boxlo,domain->boxlo,3*sizeof(double));
memcpy(prd,domain->prd,3*sizeof(double));
// store grid bounds and proc neighs
if (flag == 1) {
if (layout == REGULAR) {
// this assumes gcomm = world
int (*procneigh)[2] = comm->procneigh;
initialize(gcomm,gnx,gny,gnz,
ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
exlo,exhi,eylo,eyhi,ezlo,ezhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1],
procneigh[2][0],procneigh[2][1]);
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
exlo,exhi,eylo,eyhi,ezlo,ezhi,
procneigh[0][0],procneigh[0][1],
procneigh[1][0],procneigh[1][1],
procneigh[2][0],procneigh[2][1]);
} else {
initialize(gcomm,gnx,gny,gnz,
ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
exlo,exhi,eylo,eyhi,ezlo,ezhi,
0,0,0,0,0,0);
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
exlo,exhi,eylo,eyhi,ezlo,ezhi,
0,0,0,0,0,0);
}
} else if (flag == 2) {
if (layout == REGULAR) {
initialize(gcomm,gnx,gny,gnz,
ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
exlo,exhi,eylo,eyhi,ezlo,ezhi);
store(ixlo,ixhi,iylo,iyhi,izlo,izhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
oxlo,oxhi,oylo,oyhi,ozlo,ozhi,
exlo,exhi,eylo,eyhi,ezlo,ezhi);
} else {
error->all(FLERR,"Grid3d does not support tiled layout with neighbor procs");
}
@ -161,28 +306,18 @@ Grid3d::~Grid3d()
}
/* ----------------------------------------------------------------------
store constructor args in local variables
store grid bounds and proc neighs in local variables
------------------------------------------------------------------------- */
void Grid3d::initialize(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,
int fxlo, int fxhi, int fylo, int fyhi,
int fzlo, int fzhi,
int pxlo, int pxhi, int pylo, int pyhi,
int pzlo, int pzhi)
void Grid3d::store(int ixlo, int ixhi, int iylo, int iyhi,
int izlo, int izhi,
int oxlo, int oxhi, int oylo, int oyhi,
int ozlo, int ozhi,
int fxlo, int fxhi, int fylo, int fyhi,
int fzlo, int fzhi,
int pxlo, int pxhi, int pylo, int pyhi,
int pzlo, int pzhi)
{
gridcomm = gcomm;
MPI_Comm_rank(gridcomm,&me);
MPI_Comm_size(gridcomm,&nprocs);
nx = gnx;
ny = gny;
nz = gnz;
inxlo = ixlo;
inxhi = ixhi;
inylo = iylo;
@ -229,7 +364,7 @@ void Grid3d::initialize(MPI_Comm gcomm,
/* ---------------------------------------------------------------------- */
void Grid3d::query_global_size(int &nxgrid, int &nygrid, int &nzgrid)
void Grid3d::query_size(int &nxgrid, int &nygrid, int &nzgrid)
{
nxgrid = nx;
nygrid = ny;
@ -238,7 +373,8 @@ void Grid3d::query_global_size(int &nxgrid, int &nygrid, int &nzgrid)
/* ---------------------------------------------------------------------- */
void Grid3d::query_in_bounds(int &xlo, int &xhi, int &ylo, int &yhi, int &zlo, int &zhi)
void Grid3d::query_bounds(int &xlo, int &xhi, int &ylo, int &yhi,
int &zlo, int &zhi)
{
xlo = inxlo;
xhi = inxhi;
@ -250,6 +386,14 @@ void Grid3d::query_in_bounds(int &xlo, int &xhi, int &ylo, int &yhi, int &zlo, i
/* ---------------------------------------------------------------------- */
void Grid3d::query_box(int dim, double &lo, double &delta)
{
lo = boxlo[dim];
delta = prd[dim] / ngrid[dim];
}
/* ---------------------------------------------------------------------- */
void Grid3d::setup(int &nbuf1, int &nbuf2)
{
if (layout == REGULAR) setup_regular(nbuf1,nbuf2);

25
src/grid3d.h
View File

@ -22,13 +22,18 @@ class Grid3d : protected Pointers {
public:
enum { KSPACE = 0, PAIR = 1, FIX = 2 }; // calling classes
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(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, 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 query_global_size(int &, int &, int &);
void query_in_bounds(int &, int &, int &, int &, int &, int &);
void query_size(int &, int &, int &);
void query_bounds(int &, int &, int &, int &, int &, int &);
void query_box(int, double &, double &);
void setup(int &, int &);
int ghost_adjacent();
void forward_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
@ -41,6 +46,10 @@ class Grid3d : protected Pointers {
MPI_Comm gridcomm; // communicator for this class
// usually world, but MSM calls with subset
int ngrid[3]; // global grid size
double boxlo[3]; // current box that grid is mapped to
double prd[3];
// inputs from caller via constructor
int nx, ny, nz; // size of global grid in all 3 dims
@ -176,8 +185,8 @@ class Grid3d : protected Pointers {
// internal methods
// -------------------------------------------
void initialize(MPI_Comm, int, int, int, 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 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);
virtual void setup_regular(int &, int &);
virtual void setup_tiled(int &, int &);
void ghost_box_drop(int *, int *);