Hello All,
I use LAMMPS (SPH Package from the 8 April 2021 version) for modeling electrodeposition and I am currently attempting to implement a connectivity check of particles. My system has a type 1 and type 2 particle. As the simulation cycles some type 2 particles will break off from the origional portion of type 2 particles and become isolated (i.e. surrounded by type 1 particles). I am using a breadth-first search in order to determine the connectivity of particles and this seems to work properly within each processor domain. However, when communicating with other processors, there seems to be an issue. The code is below and the issue can be visualized in the attached image. The biggest issue is that the boundary particles are all being turned to type 3 when they shouldn’t be.
Here is the code:
/* ---------------------------------------------------------------------- */
FixSPHconnectivityBFS::FixSPHconnectivityBFS(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if ((atom->esph_flag != 1) || (atom->rho_flag != 1))
error->all(FLERR,
"fix sph/precipitation/dissolution/liquidRC/Connectivity/BFS command requires atom_style with both energy and density, e.g. meso");
if (narg != 6)
error->all(FLERR,"Illegal number of arguments for fix sph/connectivity/BFS command");
time_integrate = 0;
// Obtain the neighbor cutoff distance
int m = 3;
dx = atof(arg[m++]);
cCeq = atof(arg[m++]);
cAeq = atof(arg[m++]);
if (dx <= 0)
error->all(FLERR,"Illegal value for equilibrium cutoff distance");
if (cCeq <= 0)
error->all(FLERR,"Illegal value for equilibrium cation concentration");
if (cAeq <= 0)
error->all(FLERR,"Illegal value for equilibrium anion concentration");
// find the concentration property
int fcC;
int icC = atom->find_custom("cC", fcC);
if (icC < 0)
error->all(FLERR,
"Can't find property cC for fix sph/connectivity/BFS");
cC = atom->dvector[icC];
// find the concentration property
int fcA;
int icA = atom->find_custom("cA", fcA);
if (icA < 0)
error->all(FLERR,
"Can't find property cA for fix sph/connectivity/BFS");
cA = atom->dvector[icA];
// find the concentration property
int fdcA;
int idcA = atom->find_custom("dcA", fdcA);
if (idcA < 0)
error->all(FLERR,
"Can't find property dcA for fix sph/connectivity/BFS");
dcA = atom->dvector[idcA];
// find the mass of A property
int fmM;
int imM = atom->find_custom("mM", fmM);
if (imM < 0)
error->all(FLERR,
"Can't find property mM for fix sph/connectivity/BFS");
mM = atom->dvector[imM];
// find the potential property
int flocal_pot;
int ilocal_pot = atom->find_custom("local_pot", flocal_pot);
if (ilocal_pot < 0)
error->all(FLERR,
"Can't find property local_pot for fix sph/connectivity/BFS");
local_pot = atom->dvector[ilocal_pot];
// find the anode particles
int fisAnode;
int iisAnode = atom->find_custom("isAnode", fisAnode);
if (iisAnode < 0)
error->all(FLERR,
"Can't find property isAnode for fix sph/connectivity/BFS");
isAnode = atom->dvector[iisAnode];
// set comm size by this fix
comm_forward = 2;
comm_reverse = 2;
}
/* ---------------------------------------------------------------------- */
FixSPHconnectivityBFS::~FixSPHconnectivityBFS() {
}
/* ---------------------------------------------------------------------- */
int FixSPHconnectivityBFS::setmask() {
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSPHconnectivityBFS::init() {
// need a full neighbor list, built whenever re-neighboring occurs
int irequest = neighbor->request((void *) this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
dtcA = update->dt;
}
/* ---------------------------------------------------------------------- */
void FixSPHconnectivityBFS::init_list(int, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void FixSPHconnectivityBFS::end_of_step()
{
int i, j, ii, jj, itype, jtype, jshortest, inum, jnum;
int *ilist, *jlist, *numneigh, **firstneigh;
double delx, dely, delz;
double shortest, rsq, r, cAdist;
double tempDistance;
double **x = atom->x;
double **v = atom->v;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Create memory for connectivity array
memory->create(connected, nall, "fix:connected");
// Communicate the local atom type to the ghost atoms
comm->forward_comm_fix(this);
for (i = 0; i < nall; i++) {
connected[i] = 0.0;
}
std::queue<int> queue;
//First pass: mark all anode particles and add them to the queue
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
if(itype == 2 && isAnode[i] == 1 && i < nlocal) {
connected[i] = 1;
queue.push(i); // Add directly connected anode particle to the queue
}
}
//Use BFS to propogate connectivity
while(!queue.empty()) {
i = queue.front();
queue.pop();
if (i >= nlocal) continue;
jnum = numneigh[i];
jlist = firstneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
// Check if the distance between particles i and j is within an acceptable range
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
r = sqrt(delx*delx + dely*dely + delz*delz);
if (jtype == 2 && connected[j] == 0 && r <= 1.75*dx) {
connected[j] = 1;
queue.push(j); // Add newly connected particles to the queue
}
}
}
// Send the phase change status
comm->reverse_comm_fix(this);
for (i = 0; i < nlocal; i++) {
if (connected[i] == 0 && type[i] == 2) {
type[i] = 3;
}
}
// Destroy the memory created for ischangecC
memory->destroy(connected);
}
/* ---------------------------------------------------------------------- */
int FixSPHconnectivityBFS::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = atom->type[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSPHconnectivityBFS::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
atom->type[i] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int FixSPHconnectivityBFS::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
int *type = atom->type;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = connected[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSPHconnectivityBFS::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
int *type = atom->type;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
connected[j] += buf[m++];
}
}
Note: I have also attempted to add the atom type to the reverse communication and send the reverse communication following the phase change.
In the above figure the gray particles are type 2, the blue are type 1 and the black are type 3. The red arrows are added to show where the communication seems to be failing. They system is running on 32 processors.
Does anyone understand where I may be going wrong in my communication?
Thanks!