more on shift factors

This commit is contained in:
Steve Plimpton
2022-10-28 18:03:08 -06:00
parent 861e3b5876
commit 19fad284af
2 changed files with 16 additions and 34 deletions

View File

@ -41,13 +41,8 @@ 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 to match fix ttm
// also it allows fix ave/chunk to spatially average consistently
static constexpr int OFFSET = 16384;
static constexpr double SHIFT = 0.5;
/* ---------------------------------------------------------------------- */
@ -61,7 +56,6 @@ FixTTMGrid::FixTTMGrid(LAMMPS *lmp, int narg, char **arg) :
if (outfile) error->all(FLERR,"Fix ttm/grid does not support outfile option - "
"use dump grid command or restart files instead");
shift = OFFSET + SHIFT;
skin_original = neighbor->skin;
}
@ -142,9 +136,9 @@ void FixTTMGrid::post_force(int /*vflag*/)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET;
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + OFFSET) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + OFFSET) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + OFFSET) - OFFSET;
// flag if ix,iy,iz is not within my ghost cell range
@ -201,24 +195,18 @@ void FixTTMGrid::end_of_step()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET;
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + OFFSET) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + OFFSET) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + OFFSET) - OFFSET;
net_energy_transfer[iz][iy][ix] +=
(flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2]);
}
printf("AAA Telec000 %g Net000 %g\n",
T_electron[0][0][0],net_energy_transfer[0][0][0]);
grid->reverse_comm(Grid3d::FIX,this,1,sizeof(double),0,
grid_buf1,grid_buf2,MPI_DOUBLE);
printf("BBB Telec000 %g Net000 %g\n",
T_electron[0][0][0],net_energy_transfer[0][0][0]);
// clang-format off
// num_inner_timesteps = # of inner steps (thermal solves)
@ -272,9 +260,6 @@ void FixTTMGrid::end_of_step()
grid->forward_comm(Grid3d::FIX,this,1,sizeof(double),0,
grid_buf1,grid_buf2,MPI_DOUBLE);
}
printf("CCC Telec000 %g Net000 %g\n",
T_electron[0][0][0],net_energy_transfer[0][0][0]);
}
/* ----------------------------------------------------------------------
@ -500,7 +485,8 @@ void FixTTMGrid::reset_grid()
int tmp[12];
double maxdist = 0.5 * neighbor->skin;
Grid3d *gridnew = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, maxdist, 1, SHIFT,
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]);
@ -646,7 +632,7 @@ void FixTTMGrid::allocate_grid()
{
double maxdist = 0.5 * neighbor->skin;
grid = new Grid3d(lmp, world, nxgrid, nygrid, nzgrid, maxdist, 1, SHIFT,
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);

View File

@ -39,11 +39,8 @@ enum{ONE,RUNNING,WINDOW};
enum{DISCARD,KEEP};
// 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
static constexpr int OFFSET = 16384;
static constexpr double SHIFT = 0.0;
/* ---------------------------------------------------------------------- */
@ -398,7 +395,7 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
else if (modegrid) maxdist = 0.0;
if (dimension == 2) {
grid2d = new Grid2d(lmp, world, nxgrid, nygrid, maxdist, 0, SHIFT,
grid2d = new Grid2d(lmp, world, nxgrid, nygrid, maxdist, 0, 0.5,
nxlo_in, nxhi_in, nylo_in, nyhi_in,
nxlo_out, nxhi_out, nylo_out, nyhi_out);
@ -414,7 +411,7 @@ 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, SHIFT,
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);
@ -846,7 +843,6 @@ void FixAveGrid::atom2grid()
if (triclinic) domain->x2lamda(nlocal);
int flag = 0;
double shift = OFFSET + SHIFT;
if (dimension == 2) {
for (i = 0; i < nlocal; i++) {
@ -855,8 +851,8 @@ void FixAveGrid::atom2grid()
continue;
}
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET;
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + OFFSET) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + OFFSET) - OFFSET;
if (ix < nxlo_out || ix > nxhi_out) {
if (periodicity[0]) flag = 1;
@ -890,9 +886,9 @@ void FixAveGrid::atom2grid()
continue;
}
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + shift) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + shift) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + shift) - OFFSET;
ix = static_cast<int> ((x[i][0]-boxlo[0])*dxinv + OFFSET) - OFFSET;
iy = static_cast<int> ((x[i][1]-boxlo[1])*dyinv + OFFSET) - OFFSET;
iz = static_cast<int> ((x[i][2]-boxlo[2])*dzinv + OFFSET) - OFFSET;
if (ix < nxlo_out || ix > nxhi_out) {
if (periodicity[0]) flag = 1;