Commit 0512e788 authored by Tim Mattox's avatar Tim Mattox
Browse files

USER-DPD: new neighbor list code for SSA that gives neighbors to ghosts.

This simplifies the processing of the neighbor list in fix_shardlow.
NOTE: pair evaluation order changes, causing numerical differences!
parent 2db66e49
Loading
Loading
Loading
Loading
+11 −7
Original line number Diff line number Diff line
@@ -148,6 +148,7 @@ void FixShardlow::init()
  int irequest = neighbor->request(this,instance_me);
  neighbor->requests[irequest]->pair = 0;
  neighbor->requests[irequest]->fix  = 1;
  neighbor->requests[irequest]->ghost= 1;
  neighbor->requests[irequest]->ssa  = 1;
}

@@ -498,7 +499,7 @@ void FixShardlow::ssa_update_dpde(

void FixShardlow::initial_integrate(int vflag)
{
  int i,ii,inum;
  int i,ii,inum,anum;
  int *ilist;

  int nlocal = atom->nlocal;
@@ -531,10 +532,12 @@ void FixShardlow::initial_integrate(int vflag)
  v_t0 = (double (*)[3]) memory->smalloc(sizeof(double)*3*nghost, "FixShardlow:v_t0");

  inum = list->inum;
  anum = inum + list->gnum;
  ilist = list->ilist;

  dtsqrt = sqrt(update->dt);

  ii = 0;
  //Loop over all 14 directions (8 stages)
  for (airnum = 1; airnum <=8; airnum++){

@@ -549,15 +552,16 @@ void FixShardlow::initial_integrate(int vflag)
      }
    }

    // Loop over neighbors of my atoms
    for (ii = 0; ii < inum; ii++) {
    // process neighbors in this AIR
    while (ii < anum) {
      i = ilist[ii];
      int start = (airnum < 2) ? 0 : list->ndxAIR_ssa[i][airnum - 2];
      int len = list->ndxAIR_ssa[i][airnum - 1] - start;
      if (atom->ssaAIR[i] > airnum) break; /* done with curent AIR */
      int len = list->numneigh[i];
      if (len > 0) {
        if (useDPDE) ssa_update_dpde(i, &(list->firstneigh[i][start]), len);
        else ssa_update_dpd(i, &(list->firstneigh[i][start]), len);
        if (useDPDE) ssa_update_dpde(i, &(list->firstneigh[i][0]), len);
        else ssa_update_dpd(i, &(list->firstneigh[i][0]), len);
      }
      ii++;
    }

    // Communicate the ghost deltas to the atom owners
+13 −14
Original line number Diff line number Diff line
@@ -33,14 +33,13 @@ NBinSSA::NBinSSA(LAMMPS *lmp) : NBinStandard(lmp)
  bins_ssa = NULL;
  maxhead_ssa = 0;
  binhead_ssa = NULL;
  gbinhead_ssa = NULL;
  for (int i = 0; i < 9; i++) gairhead_ssa[i] = -1;
}

NBinSSA::~NBinSSA()
{
  memory->destroy(bins_ssa);
  memory->destroy(binhead_ssa);
  memory->destroy(gbinhead_ssa);
}

/* ----------------------------------------------------------------------
@@ -62,8 +61,11 @@ void NBinSSA::bin_atoms()

  last_bin = update->ntimestep;

  for (i = 0; i < 9; i++) {
    gairhead_ssa[i] = -1;
  }

  for (i = 0; i < mbins; i++) {
    gbinhead_ssa[i] = -1;
    binhead_ssa[i] = -1;
  }

@@ -73,19 +75,19 @@ void NBinSSA::bin_atoms()
    int bitmask = group->bitmask[includegroup];
    int nowned = atom->nlocal; // NOTE: nlocal was set to atom->nfirst above
    for (i = nall-1; i >= nowned; i--) {
      if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
      ibin = ssaAIR[i];
      if (ibin < 2) continue; // skip ghost atoms not in AIR
      if (mask[i] & bitmask) {
        ibin = coord2bin(x[i]);
        bins_ssa[i] = gbinhead_ssa[ibin];
        gbinhead_ssa[ibin] = i;
        bins_ssa[i] = gairhead_ssa[ibin];
        gairhead_ssa[ibin] = i;
      }
    }
  } else {
    for (i = nall-1; i >= nlocal; i--) {
      if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
      ibin = coord2bin(x[i]);
      bins_ssa[i] = gbinhead_ssa[ibin];
      gbinhead_ssa[ibin] = i;
      ibin = ssaAIR[i];
      if (ibin < 2) continue; // skip ghost atoms not in AIR
      bins_ssa[i] = gairhead_ssa[ibin];
      gairhead_ssa[ibin] = i;
    }
  }
  for (i = nlocal-1; i >= 0; i--) {
@@ -103,10 +105,8 @@ void NBinSSA::bin_atoms_setup(int nall)

  if (mbins > maxhead_ssa) {
    maxhead_ssa = mbins;
    memory->destroy(gbinhead_ssa);
    memory->destroy(binhead_ssa);
    memory->create(binhead_ssa,maxhead_ssa,"binhead_ssa");
    memory->create(gbinhead_ssa,maxhead_ssa,"gbinhead_ssa");
  }

  if (nall > maxbin_ssa) {
@@ -125,7 +125,6 @@ bigint NBinSSA::memory_usage()
  if (maxbin_ssa) bytes += memory->usage(bins_ssa,maxbin_ssa);
  if (maxhead_ssa) {
    bytes += memory->usage(binhead_ssa,maxhead_ssa);
    bytes += memory->usage(gbinhead_ssa,maxhead_ssa);
  }
  return bytes;
}
+1 −1
Original line number Diff line number Diff line
@@ -32,7 +32,7 @@ class NBinSSA : public NBinStandard {
  int *bins_ssa;             // index of next atom in each bin
  int maxbin_ssa;            // size of bins_ssa array
  int *binhead_ssa;          // index of 1st local atom in each bin
  int *gbinhead_ssa;         // index of 1st ghost atom in each bin
  int gairhead_ssa[9];       // index of 1st ghost atom in each AIR
  int maxhead_ssa;           // size of binhead_ssa and gbinhead_ssa arrays

  NBinSSA(class LAMMPS *);
+68 −79
Original line number Diff line number Diff line
@@ -32,12 +32,6 @@

using namespace LAMMPS_NS;

// allocate space for static class variable
// prototype for non-class function

static int *ssaAIRptr;
static int cmp_ssaAIR(const void *, const void *);

/* ---------------------------------------------------------------------- */

NPairHalfBinNewtonSSA::NPairHalfBinNewtonSSA(LAMMPS *lmp) : NPair(lmp) {}
@@ -64,9 +58,7 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
  tagint **special = atom->special;
  int **nspecial = atom->nspecial;
  int nlocal = atom->nlocal;
  int nall = nlocal + atom->nghost;
  if (includegroup) nlocal = atom->nfirst;
  int *ssaAIR = atom->ssaAIR;

  int *molindex = atom->molindex;
  int *molatom = atom->molatom;
@@ -89,16 +81,18 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
  if (!nb_ssa) error->one(FLERR, "NBin wasn't a NBinSSA object");
  int *bins_ssa = nb_ssa->bins_ssa;
  int *binhead_ssa = nb_ssa->binhead_ssa;
  int *gbinhead_ssa = nb_ssa->gbinhead_ssa;
  int *gairhead_ssa = &(nb_ssa->gairhead_ssa[0]);

  int inum = 0;
  int gnum = 0;
  int xbin,ybin,zbin,xbin2,ybin2,zbin2;
  int **stencilxyz = ns_ssa->stencilxyz;

  ipage->reset();

  // loop over owned atoms, storing half of the neighbors

  for (i = 0; i < nlocal; i++) {
    int AIRct[8] = { 0 };
    n = 0;
    neighptr = ipage->vget();

@@ -175,16 +169,47 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
        }
      }
    }
    AIRct[0] = n;

    ilist[inum++] = i;
    firstneigh[i] = neighptr;
    numneigh[i] = n;
    ipage->vgot(n);
    if (ipage->status())
      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
  }

  list->inum = inum;

  // loop over AIR ghost atoms, storing their local neighbors
  // since these are ghosts, must check if stencil bin is out of bounds
  for (int airnum = 2; airnum <= 8; airnum++) {
    for (i = gairhead_ssa[airnum]; i >= 0; i = bins_ssa[i]) {
      n = 0;
      neighptr = ipage->vget();

      itype = type[i];
      xtmp = x[i][0];
      ytmp = x[i][1];
      ztmp = x[i][2];
      if (moltemplate) {
        imol = molindex[i];
        iatom = molatom[i];
        tagprev = tag[i] - iatom - 1;
      }

      ibin = coord2bin(x[i],xbin,ybin,zbin);

      // loop over AIR ghost atoms in all bins in "full" stencil
      // Note: the non-AIR ghost atoms have already been filtered out
    // That is a significant time savings because of the "full" stencil
    // Note2: only non-pure locals can have ghosts as neighbors

    if (ssaAIR[i] == 1) for (k = 0; k < nstencil_full; k++) {
      for (j = gbinhead_ssa[ibin+stencil[k]]; j >= 0;
           j = bins_ssa[j]) {
      for (k = 0; k < nstencil_full; k++) {
        xbin2 = xbin + stencilxyz[k][0];
        ybin2 = ybin + stencilxyz[k][1];
        zbin2 = zbin + stencilxyz[k][2];
        // since we only care about ghost to local neighbors, these "bounds" could be inset
        if (xbin2 < 0 || xbin2 >= mbinx ||
            ybin2 < 0 || ybin2 >= mbiny ||
            zbin2 < 0 || zbin2 >= mbinz) continue;
        for (j = binhead_ssa[ibin+stencil[k]]; j >= 0; j = bins_ssa[j]) {

          jtype = type[j];
          if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
@@ -203,58 +228,22 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
                                     onemols[imol]->nspecial[iatom],
                                     tag[j]-tagprev);
              else which = 0;
            if (which == 0) {
              neighptr[n++] = j;
              ++(AIRct[ssaAIR[j] - 1]);
            } else if (domain->minimum_image_check(delx,dely,delz)) {
              neighptr[n++] = j;
              ++(AIRct[ssaAIR[j] - 1]);
            } else if (which > 0) {
              neighptr[n++] = j ^ (which << SBBITS);
              ++(AIRct[ssaAIR[j] - 1]);
            }
          } else {
              if (which == 0) neighptr[n++] = j;
              else if (domain->minimum_image_check(delx,dely,delz))
                neighptr[n++] = j;
            ++(AIRct[ssaAIR[j] - 1]);
          }
              else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
            } else neighptr[n++] = j;
          }
        }
      }

    ilist[inum++] = i;
      if (n > 0) ilist[inum + (gnum++)] = i;
      firstneigh[i] = neighptr;
      numneigh[i] = n;
      ipage->vgot(n);
      if (ipage->status())
      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");

    // sort the ghosts in the neighbor list by their ssaAIR number

    ssaAIRptr = atom->ssaAIR;
    qsort(&(neighptr[AIRct[0]]), n - AIRct[0], sizeof(int), cmp_ssaAIR);

    // do a prefix sum on the counts to turn them into indexes

    list->ndxAIR_ssa[i][0] = AIRct[0];
    for (int ndx = 1; ndx < 8; ++ndx) {
      list->ndxAIR_ssa[i][ndx] = AIRct[ndx] + list->ndxAIR_ssa[i][ndx - 1];
        error->one(FLERR,"Neighbor (ghost) list overflow, boost neigh_modify one");
    }
  }

  list->inum = inum;
}

/* ----------------------------------------------------------------------
   comparison function invoked by qsort()
   accesses static class member ssaAIRptr, set before call to qsort()
------------------------------------------------------------------------- */

static int cmp_ssaAIR(const void *iptr, const void *jptr)
{
  int i = NEIGHMASK & *((int *) iptr);
  int j = NEIGHMASK & *((int *) jptr);
  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;
  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;
  return 0;
  list->gnum = gnum;
}
+1 −1
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@

NPairStyle(half/bin/newton/ssa,
           NPairHalfBinNewtonSSA,
           NP_HALF | NP_BIN | NP_NEWTON | NP_ORTHO | NP_SSA)
           NP_HALF | NP_BIN | NP_NEWTON | NP_ORTHO | NP_SSA | NP_GHOST)

#else

Loading