more debugging of restarts

This commit is contained in:
Steve Plimpton
2021-08-25 16:28:40 -06:00
parent ae94a60d4a
commit a0dfae9876
2 changed files with 41 additions and 40 deletions

View File

@ -287,16 +287,18 @@ void FixTTM::post_force(int /*vflag*/)
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *boxlo = domain->boxlo;
double dxinv = nxgrid/domain->xprd;
double dyinv = nygrid/domain->yprd;
double dzinv = nzgrid/domain->zprd;
// apply damping and thermostat to all atoms in fix group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
ix = static_cast<int>(xscale*nxgrid + shift) - OFFSET;
iy = static_cast<int>(yscale*nygrid + shift) - OFFSET;
iz = static_cast<int>(zscale*nzgrid + shift) - OFFSET;
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;
if (ix < 0) ix += nxgrid;
if (iy < 0) iy += nygrid;
if (iz < 0) iz += nzgrid;
@ -354,6 +356,11 @@ void FixTTM::end_of_step()
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *boxlo = domain->boxlo;
double dxinv = nxgrid/domain->xprd;
double dyinv = nygrid/domain->yprd;
double dzinv = nzgrid/domain->zprd;
for (iz = 0; iz < nzgrid; iz++)
for (iy = 0; iy < nygrid; iy++)
for (ix = 0; ix < nxgrid; ix++)
@ -361,25 +368,19 @@ void FixTTM::end_of_step()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
ix = static_cast<int>(xscale*nxgrid + shift) - OFFSET;
iy = static_cast<int>(yscale*nygrid + shift) - OFFSET;
iz = static_cast<int>(zscale*nzgrid + shift) - OFFSET;
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;
if (ix < 0) ix += nxgrid;
if (iy < 0) iy += nygrid;
if (iz < 0) iz += nzgrid;
if (ix >= nxgrid) ix -= nxgrid;
if (iy >= nygrid) iy -= nygrid;
if (iz >= nzgrid) iz -= nzgrid;
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("FLANG i %d xyz %g %g %g\n",i,
// flangevin[i][0],
// flangevin[i][1],
// flangevin[i][2]);
}
outflag = 0;
@ -443,15 +444,12 @@ void FixTTM::end_of_step()
inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity *
((T_electron_old[iz][iy][xright] +
T_electron_old[iz][iy][xleft] -
2*T_electron_old[iz][iy][ix])/dx/dx +
(T_electron_old[iz][yright][ix] +
T_electron_old[iz][yleft][ix] -
2*T_electron_old[iz][iy][ix])/dy/dy +
(T_electron_old[zright][iy][ix] +
T_electron_old[zleft][iy][ix] -
2*T_electron_old[iz][iy][ix])/dz/dz) -
((T_electron_old[iz][iy][xright] + T_electron_old[iz][iy][xleft] -
2.0*T_electron_old[iz][iy][ix])/dx/dx +
(T_electron_old[iz][yright][ix] + T_electron_old[iz][yleft][ix] -
2.0*T_electron_old[iz][iy][ix])/dy/dy +
(T_electron_old[zright][iy][ix] + T_electron_old[zleft][iy][ix] -
2.0*T_electron_old[iz][iy][ix])/dz/dz) -
(net_energy_transfer_all[iz][iy][ix])/del_vol);
}

View File

@ -79,8 +79,13 @@ void FixTTMGrid::post_constructor()
ngridout*sizeof(double));
// set initial electron temperatures from user input file
// communicate new T_electron values to ghost grid points
if (infile) read_electron_temperatures(infile);
if (infile) {
read_electron_temperatures(infile);
gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0,
gc_buf1,gc_buf2,MPI_DOUBLE);
}
}
/* ---------------------------------------------------------------------- */
@ -118,7 +123,6 @@ 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;
@ -183,6 +187,7 @@ void FixTTMGrid::end_of_step()
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;
net_energy_transfer[iz][iy][ix] +=
(flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2]);
@ -228,21 +233,21 @@ void FixTTMGrid::end_of_step()
(electronic_thermal_conductivity *
((T_electron_old[iz][iy][ix-1] + T_electron_old[iz][iy][ix+1] -
2*T_electron_old[iz][iy][ix])*dxinv*dxinv +
2.0*T_electron_old[iz][iy][ix])*dxinv*dxinv +
(T_electron_old[iz][iy-1][ix] + T_electron_old[iz][iy+1][ix] -
2*T_electron_old[iz][iy][ix])*dyinv*dyinv +
2.0*T_electron_old[iz][iy][ix])*dyinv*dyinv +
(T_electron_old[iz-1][iy][ix] + T_electron_old[iz+1][iy][ix] -
2*T_electron_old[iz][iy][ix])*dzinv*dzinv) -
2.0*T_electron_old[iz][iy][ix])*dzinv*dzinv) -
net_energy_transfer[iz][iy][ix]/volgrid);
}
// communicate new T_electron values to ghost grid points
gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0,
gc_buf1,gc_buf2,MPI_DOUBLE);
}
// communicate new T_electron values to ghost grid points
gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0,
gc_buf1,gc_buf2,MPI_DOUBLE);
// output of grid temperatures to file
if (outfile && (update->ntimestep % outevery == 0)) {
@ -364,11 +369,6 @@ void FixTTMGrid::read_electron_temperatures(const char *filename)
error->all(FLERR,"Fix ttm/grid infile did not set all temperatures");
memory->destroy3d_offset(T_initial_set,nzlo_in,nylo_in,nxlo_in);
// communicate new T_electron values to ghost grid points
gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0,
gc_buf1,gc_buf2,MPI_DOUBLE);
}
/* ----------------------------------------------------------------------
@ -593,6 +593,9 @@ void FixTTMGrid::restart(char *buf)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
iglobal = nygrid*nxgrid*iz + nxgrid*iy + ix;
T_electron[iz][iy][ix] = rlist[n+iglobal];
if (ix == 0 && iy == 0 & iz == 0)
printf("READRESTART 000 n+iglobal %d %d value %g\n",n,iglobal,
T_electron[ix][iy][iz]);
}
}