Commit 5fb5f70e authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

update USER-OMP code for shear history neighbor refactoring

parent 53aa92cf

src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp

+9 −68
Original line number Diff line number Diff line
@@ -18,9 +18,6 @@ 
#include "neigh_list.h"

#include "atom.h"

#include "atom_vec.h"

#include "molecule.h"

#include "domain.h"

#include "fix_shear_history.h"

#include "my_page.h"

#include "error.h"
@@ -34,7 +31,6 @@ NPairHalfSizeBinNewtoffOmp::NPairHalfSizeBinNewtoffOmp(LAMMPS *lmp) : 
/* ----------------------------------------------------------------------

   size particles

   binned neighbor list construction with partial Newton's 3rd law

   shear history must be accounted for when a neighbor pair is added

   each owned atom i checks own bin and surrounding bins in non-Newton stencil

   pair stored once if i,j are both owned and i < j

   pair stored by me if j is ghost (also stored by proc owning j)
@@ -43,30 +39,20 @@ NPairHalfSizeBinNewtoffOmp::NPairHalfSizeBinNewtoffOmp(LAMMPS *lmp) : 
void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)

{

  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;



  FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history;

  NeighList * listhistory = list->listhistory;

  const int history = list->history;

  const int mask_history = 3 << SBBITS;



  NPAIR_OMP_INIT;



#if defined(_OPENMP)

#pragma omp parallel default(none) shared(list,listhistory)

#pragma omp parallel default(none) shared(list)

#endif

  NPAIR_OMP_SETUP(nlocal);



  int i,j,k,m,n,nn,ibin,dnum,dnumbytes;

  int i,j,k,m,n,nn,ibin;

  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;

  double radi,radsum,cutsq;

  int *neighptr,*touchptr;

  double *shearptr;

  MyPage<int> *ipage_touch;

  MyPage<double> *dpage_shear;



  int *npartner;

  tagint **partner;

  double **shearpartner;

  int **firsttouch;

  double **firstshear;

  int *neighptr;



  // loop over each atom, storing neighbors
@@ -85,29 +71,10 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) 
  MyPage<int> &ipage = list->ipage[tid];

  ipage.reset();



  if (fix_history) {

    npartner = fix_history->npartner;

    partner = fix_history->partner;

    shearpartner = fix_history->shearpartner;

    firsttouch = listhistory->firstneigh;

    firstshear = listhistory->firstdouble;

    ipage_touch = listhistory->ipage+tid;

    dpage_shear = listhistory->dpage+tid;

    dnum = listhistory->dnum;

    dnumbytes = dnum * sizeof(double);

    ipage_touch->reset();

    dpage_shear->reset();

  }



  for (i = ifrom; i < ito; i++) {



    n = 0;

    neighptr = ipage.vget();

    if (fix_history) {

      nn = 0;

      touchptr = ipage_touch->vget();

      shearptr = dpage_shear->vget();

    }



    xtmp = x[i][0];

    ytmp = x[i][1];
@@ -133,29 +100,10 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) 
        cutsq = (radsum+skin) * (radsum+skin);



        if (rsq <= cutsq) {

          neighptr[n] = j;



          if (fix_history) {

            if (rsq < radsum*radsum) {

              for (m = 0; m < npartner[i]; m++)

                if (partner[i][m] == tag[j]) break;

              if (m < npartner[i]) {

                touchptr[n] = 1;

                memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);

                nn += dnum;

              } else {

                touchptr[n] = 0;

                memcpy(&shearptr[nn],zeroes,dnumbytes);

                nn += dnum;

              }

            } else {

              touchptr[n] = 0;

              memcpy(&shearptr[nn],zeroes,dnumbytes);

              nn += dnum;

            }

          }



          n++;

          if (history && rsq < radsum*radsum)

            neighptr[n++] = j ^ mask_history;

          else

            neighptr[n++] = j;

        }

      }

    }
@@ -166,13 +114,6 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) 
    ipage.vgot(n);

    if (ipage.status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");



    if (fix_history) {

      firsttouch[i] = touchptr;

      firstshear[i] = shearptr;

      ipage_touch->vgot(n);

      dpage_shear->vgot(nn);

    }

  }

  NPAIR_OMP_CLOSE;

  list->inum = nlocal;

src/USER-OMP/npair_half_size_bin_newton_omp.cpp

+13 −97
Original line number Diff line number Diff line
@@ -18,9 +18,6 @@ 
#include "neigh_list.h"

#include "atom.h"

#include "atom_vec.h"

#include "molecule.h"

#include "domain.h"

#include "fix_shear_history.h"

#include "my_page.h"

#include "error.h"
@@ -34,7 +31,6 @@ NPairHalfSizeBinNewtonOmp::NPairHalfSizeBinNewtonOmp(LAMMPS *lmp) : 
/* ----------------------------------------------------------------------

   size particles

   binned neighbor list construction with full Newton's 3rd law

   shear history must be accounted for when a neighbor pair is added

   each owned atom i checks its own bin and other bins in Newton stencil

   every pair stored exactly once by some processor

------------------------------------------------------------------------- */
@@ -42,36 +38,20 @@ NPairHalfSizeBinNewtonOmp::NPairHalfSizeBinNewtonOmp(LAMMPS *lmp) : 
void NPairHalfSizeBinNewtonOmp::build(NeighList *list)

{

  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;



  FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history;

  NeighList * listhistory = list->listhistory;

  if (fix_history) {

    fix_history->nlocal_neigh = nlocal;

    fix_history->nall_neigh = nlocal + atom->nghost;

  }

  const int history = list->history;

  const int mask_history = 3 << SBBITS;



  NPAIR_OMP_INIT;



#if defined(_OPENMP)

#pragma omp parallel default(none) shared(list,listhistory)

#pragma omp parallel default(none) shared(list)

#endif

  NPAIR_OMP_SETUP(nlocal);



  int i,j,k,m,n,nn,ibin,dnum,dnumbytes;

  int i,j,k,m,n,nn,ibin;

  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;

  double radi,radsum,cutsq;

  int *neighptr,*touchptr;

  double *shearptr;

  MyPage<int> *ipage_touch;

  MyPage<double> *dpage_shear;



  int *npartner;

  tagint **partner;

  double **shearpartner;

  int **firsttouch;

  double **firstshear;



  // loop over each atom, storing neighbors

  int *neighptr;



  double **x = atom->x;

  double *radius = atom->radius;
@@ -88,29 +68,10 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) 
  MyPage<int> &ipage = list->ipage[tid];

  ipage.reset();



  if (fix_history) {

    npartner = fix_history->npartner;

    partner = fix_history->partner;

    shearpartner = fix_history->shearpartner;

    firsttouch = listhistory->firstneigh;

    firstshear = listhistory->firstdouble;

    ipage_touch = listhistory->ipage+tid;

    dpage_shear = listhistory->dpage+tid;

    dnum = listhistory->dnum;

    dnumbytes = dnum * sizeof(double);

    ipage_touch->reset();

    dpage_shear->reset();

  }



  for (i = ifrom; i < ito; i++) {



    n = 0;

    neighptr = ipage.vget();

    if (fix_history) {

      nn = 0;

      touchptr = ipage_touch->vget();

      shearptr = dpage_shear->vget();

    }



    xtmp = x[i][0];

    ytmp = x[i][1];
@@ -140,29 +101,10 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) 
      cutsq = (radsum+skin) * (radsum+skin);



      if (rsq <= cutsq) {

        neighptr[n] = j;



        if (fix_history) {

          if (rsq < radsum*radsum) {

            for (m = 0; m < npartner[i]; m++)

              if (partner[i][m] == tag[j]) break;

            if (m < npartner[i]) {

              touchptr[n] = 1;

              memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);

              nn += dnum;

            } else {

              touchptr[n] = 0;

              memcpy(&shearptr[nn],zeroes,dnumbytes);

              nn += dnum;

            }

          } else {

            touchptr[n] = 0;

            memcpy(&shearptr[nn],zeroes,dnumbytes);

            nn += dnum;

          }

        }



        n++;

        if (history && rsq < radsum*radsum)

          neighptr[n++] = j ^ mask_history;

        else

          neighptr[n++] = j;

      }

    }
@@ -181,29 +123,10 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) 
        cutsq = (radsum+skin) * (radsum+skin);



        if (rsq <= cutsq) {

          neighptr[n] = j;



          if (fix_history) {

            if (rsq < radsum*radsum) {

              for (m = 0; m < npartner[i]; m++)

                if (partner[i][m] == tag[j]) break;

              if (m < npartner[i]) {

                touchptr[n] = 1;

                memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);

                nn += dnum;

              } else {

                touchptr[n] = 0;

                memcpy(&shearptr[nn],zeroes,dnumbytes);

                nn += dnum;

              }

            } else {

              touchptr[n] = 0;

              memcpy(&shearptr[nn],zeroes,dnumbytes);

              nn += dnum;

            }

          }



          n++;

          if (history && rsq < radsum*radsum)

            neighptr[n++] = j ^ mask_history;

          else

            neighptr[n++] = j;

        }

      }

    }
@@ -214,13 +137,6 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) 
    ipage.vgot(n);

    if (ipage.status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");



    if (fix_history) {

      firsttouch[i] = touchptr;

      firstshear[i] = shearptr;

      ipage_touch->vgot(n);

      dpage_shear->vgot(nn);

    }

  }

  NPAIR_OMP_CLOSE;

  list->inum = nlocal;

src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp

+8 −4
Original line number Diff line number Diff line
@@ -17,8 +17,6 @@ 
#include "neigh_list.h"

#include "atom.h"

#include "atom_vec.h"

#include "molecule.h"

#include "domain.h"

#include "my_page.h"

#include "error.h"
@@ -32,7 +30,6 @@ NPairHalfSizeBinNewtonTriOmp::NPairHalfSizeBinNewtonTriOmp(LAMMPS *lmp) : 
/* ----------------------------------------------------------------------

   size particles

   binned neighbor list construction with Newton's 3rd law for triclinic

   no shear history is allowed for this option

   each owned atom i checks its own bin and other bins in triclinic stencil

   every pair stored exactly once by some processor

------------------------------------------------------------------------- */
@@ -40,6 +37,8 @@ NPairHalfSizeBinNewtonTriOmp::NPairHalfSizeBinNewtonTriOmp(LAMMPS *lmp) : 
void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)

{

  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;

  const int history = list->history;

  const int mask_history = 3 << SBBITS;



  NPAIR_OMP_INIT;

#if defined(_OPENMP)
@@ -105,7 +104,12 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) 
        radsum = radi + radius[j];

        cutsq = (radsum+skin) * (radsum+skin);



        if (rsq <= cutsq) neighptr[n++] = j;

        if (rsq <= cutsq) {

          if (history && rsq < radsum*radsum)

            neighptr[n++] = j ^ mask_history;

          else

            neighptr[n++] = j;

        }

      }

    }

src/USER-OMP/npair_half_size_nsq_newtoff_omp.cpp

+9 −71
Original line number Diff line number Diff line
@@ -19,9 +19,6 @@ 
#include "atom.h"

#include "atom_vec.h"

#include "group.h"

#include "molecule.h"

#include "domain.h"

#include "fix_shear_history.h"

#include "my_page.h"

#include "error.h"
@@ -44,34 +41,20 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) 
{

  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;

  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;



  FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history;

  NeighList * listhistory = list->listhistory;

  if (fix_history) {

    fix_history->nlocal_neigh = nlocal;

    fix_history->nall_neigh = nlocal + atom->nghost;

  }

  const int history = list->history;

  const int mask_history = 3 << SBBITS;



  NPAIR_OMP_INIT;



#if defined(_OPENMP)

#pragma omp parallel default(none) shared(list,listhistory)

#pragma omp parallel default(none) shared(list)

#endif

  NPAIR_OMP_SETUP(nlocal);



  int i,j,m,n,nn,dnum,dnumbytes;

  int i,j,m,n,nn;

  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;

  double radi,radsum,cutsq;

  int *neighptr,*touchptr;

  double *shearptr;



  int *npartner;

  tagint **partner;

  double **shearpartner;

  int **firsttouch;

  double **firstshear;

  MyPage<int> *ipage_touch;

  MyPage<double> *dpage_shear;

  int *neighptr;



  double **x = atom->x;

  double *radius = atom->radius;
@@ -89,29 +72,10 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) 
  MyPage<int> &ipage = list->ipage[tid];

  ipage.reset();



  if (fix_history) {

    npartner = fix_history->npartner;

    partner = fix_history->partner;

    shearpartner = fix_history->shearpartner;

    firsttouch = listhistory->firstneigh;

    firstshear = listhistory->firstdouble;

    ipage_touch = listhistory->ipage+tid;

    dpage_shear = listhistory->dpage+tid;

    dnum = listhistory->dnum;

    dnumbytes = dnum * sizeof(double);

    ipage_touch->reset();

    dpage_shear->reset();

  }



  for (i = ifrom; i < ito; i++) {



    n = 0;

    neighptr = ipage.vget();

    if (fix_history) {

      nn = 0;

      touchptr = ipage_touch->vget();

      shearptr = dpage_shear->vget();

    }



    xtmp = x[i][0];

    ytmp = x[i][1];
@@ -132,29 +96,10 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) 
      cutsq = (radsum+skin) * (radsum+skin);



      if (rsq <= cutsq) {

        neighptr[n] = j;



        if (fix_history) {

          if (rsq < radsum*radsum) {

            for (m = 0; m < npartner[i]; m++)

              if (partner[i][m] == tag[j]) break;

            if (m < npartner[i]) {

              touchptr[n] = 1;

              memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);

              nn += dnum;

            } else {

              touchptr[n] = 0;

              memcpy(&shearptr[nn],zeroes,dnumbytes);

              nn += dnum;

            }

          } else {

            touchptr[n] = 0;

            memcpy(&shearptr[nn],zeroes,dnumbytes);

            nn += dnum;

          }

        }



        n++;

        if (history && rsq < radsum*radsum)

          neighptr[n++] = j ^ mask_history;

        else

          neighptr[n++] = j;

      }

    }
@@ -164,13 +109,6 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) 
    ipage.vgot(n);

    if (ipage.status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");



    if (fix_history) {

      firsttouch[i] = touchptr;

      firstshear[i] = shearptr;

      ipage_touch->vgot(n);

      dpage_shear->vgot(nn);

    }

  }

  NPAIR_OMP_CLOSE;

  list->inum = nlocal;

src/USER-OMP/npair_half_size_nsq_newton_omp.cpp

+9 −70
Original line number Diff line number Diff line
@@ -19,9 +19,6 @@ 
#include "atom.h"

#include "atom_vec.h"

#include "group.h"

#include "molecule.h"

#include "domain.h"

#include "fix_shear_history.h"

#include "my_page.h"

#include "error.h"
@@ -45,34 +42,20 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) 
{

  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;

  const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;;



  FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history;

  NeighList * listhistory = list->listhistory;

  if (fix_history) {

    fix_history->nlocal_neigh = nlocal;

    fix_history->nall_neigh = nlocal+atom->nghost;

  }

  const int history = list->history;

  const int mask_history = 3 << SBBITS;



  NPAIR_OMP_INIT;



#if defined(_OPENMP)

#pragma omp parallel default(none) shared(list,listhistory)

#pragma omp parallel default(none) shared(list)

#endif

  NPAIR_OMP_SETUP(nlocal);



  int i,j,m,n,nn,itag,jtag,dnum,dnumbytes;

  int i,j,m,n,nn,itag,jtag;

  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;

  double radi,radsum,cutsq;

  int *neighptr,*touchptr;

  double *shearptr;



  int *npartner;

  tagint **partner;

  double **shearpartner;

  int **firsttouch;

  double **firstshear;

  MyPage<int> *ipage_touch;

  MyPage<double> *dpage_shear;

  int *neighptr;



  double **x = atom->x;

  double *radius = atom->radius;
@@ -90,29 +73,10 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) 
  MyPage<int> &ipage = list->ipage[tid];

  ipage.reset();



  if (fix_history) {

    npartner = fix_history->npartner;

    partner = fix_history->partner;

    shearpartner = fix_history->shearpartner;

    firsttouch = listhistory->firstneigh;

    firstshear = listhistory->firstdouble;

    ipage_touch = listhistory->ipage+tid;

    dpage_shear = listhistory->dpage+tid;

    dnum = listhistory->dnum;

    dnumbytes = dnum * sizeof(double);

    ipage_touch->reset();

    dpage_shear->reset();

  }



  for (i = ifrom; i < ito; i++) {



    n = 0;

    neighptr = ipage.vget();

    if (fix_history) {

      nn = 0;

      touchptr = ipage_touch->vget();

      shearptr = dpage_shear->vget();

    }



    itag = tag[i];

    xtmp = x[i][0];
@@ -150,29 +114,10 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) 
      cutsq = (radsum+skin) * (radsum+skin);



      if (rsq <= cutsq) {

        neighptr[n] = j;



        if (fix_history) {

          if (rsq < radsum*radsum) {

            for (m = 0; m < npartner[i]; m++)

              if (partner[i][m] == tag[j]) break;

            if (m < npartner[i]) {

              touchptr[n] = 1;

              memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes);

              nn += dnum;

            } else {

              touchptr[n] = 0;

              memcpy(&shearptr[nn],zeroes,dnumbytes);

              nn += dnum;

            }

          } else {

            touchptr[n] = 0;

            memcpy(&shearptr[nn],zeroes,dnumbytes);

            nn += dnum;

          }

        }



        n++;

        if (history && rsq < radsum*radsum)

          neighptr[n++] = j ^ mask_history;

        else

          neighptr[n++] = j;

      }

    }
@@ -183,12 +128,6 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) 
    if (ipage.status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");



    if (fix_history) {

      firsttouch[i] = touchptr;

      firstshear[i] = shearptr;

      ipage_touch->vgot(n);

      dpage_shear->vgot(nn);

    }

  }

  NPAIR_OMP_CLOSE;

  list->inum = nlocal;

CRA Git | Maintained and supported by SUSTech CRA and CCSE