lammps/src/USER-OMP/fix_shear_history_omp.cpp

/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "string.h"
#include "stdio.h"
#include "fix_shear_history_omp.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "force.h"
#include "pair.h"
#include "update.h"
#include "modify.h"
#include "error.h"

#if defined(_OPENMP)
#include <omp.h>
#endif

using namespace LAMMPS_NS;
using namespace FixConst;

/* ----------------------------------------------------------------------
   copy shear partner info from neighbor lists to per-atom arrays
   so it can be exchanged with those atoms
------------------------------------------------------------------------- */

void FixShearHistoryOMP::pre_exchange()
{

  const int nthreads = comm->nthreads;
  maxtouch = 0;

#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
  {

#if defined(_OPENMP)
    const int tid = omp_get_thread_num();
#else
    const int tid = 0;
#endif

    int i,j,ii,jj,m,n,inum,jnum;
    int *ilist,*jlist,*numneigh,**firstneigh;
    int *touch,**firsttouch;
    double *shear,*allshear,**firstshear;

    MyPage <tagint> &ipg = ipage[tid];
    MyPage <double[3]> &dpg = dpage[tid];
    ipg.reset();
    dpg.reset();

    // 1st loop over neighbor list
    // calculate nparter for each owned atom

    tagint *tag = atom->tag;

    NeighList *list = pair->list;
    inum = list->inum;
    ilist = list->ilist;
    numneigh = list->numneigh;
    firstneigh = list->firstneigh;
    firsttouch = list->listgranhistory->firstneigh;
    firstshear = list->listgranhistory->firstdouble;

    int nlocal_neigh = 0;
    if (inum) nlocal_neigh = ilist[inum-1] + 1;

    // each thread works on a fixed chunk of local and ghost atoms.
    const int ldelta = 1 + nlocal_neigh/nthreads;
    const int lfrom = tid*ldelta;
    const int lmax = lfrom +ldelta;
    const int lto = (lmax > nlocal_neigh) ? nlocal_neigh : lmax;

    // zero npartners for all current atoms and
    // clear page data structures for this thread

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


    for (ii = 0; ii < inum; ii++) {
      i = ilist[ii];
      jlist = firstneigh[i];
      jnum = numneigh[i];
      touch = firsttouch[i];

      for (jj = 0; jj < jnum; jj++) {
        if (touch[jj]) {
          if ((i >= lfrom) && (i < lto))
            npartner[i]++;

          j = jlist[jj];
          j &= NEIGHMASK;
          if ((j >= lfrom) && (j < lto))
            npartner[j]++;
        }
      }
    }

    // get page chunks to store atom IDs and shear history for my atoms

    for (ii = 0; ii < inum; ii++) {
      i = ilist[ii];
      if ((i >= lfrom) && (i < lto)) {
        n = npartner[i];
        partner[i] = ipg.get(n);
        shearpartner[i] = dpg.get(n);
        if (partner[i] == NULL || shearpartner[i] == NULL)
          error->one(FLERR,"Shear history overflow, boost neigh_modify one");
      }
    }

    // 2nd loop over neighbor list
    // store atom IDs and shear history for my atoms
    // re-zero npartner to use as counter for all my atoms

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

    for (ii = 0; ii < inum; ii++) {
      i = ilist[ii];
      jlist = firstneigh[i];
      allshear = firstshear[i];
      jnum = numneigh[i];
      touch = firsttouch[i];

      for (jj = 0; jj < jnum; jj++) {
        if (touch[jj]) {
          shear = &allshear[3*jj];
          j = jlist[jj];
          j &= NEIGHMASK;

          if ((i >= lfrom) && (i < lto)) {
            m = npartner[i];
            partner[i][m] = tag[j];
            shearpartner[i][m][0] = shear[0];
            shearpartner[i][m][1] = shear[1];
            shearpartner[i][m][2] = shear[2];
            npartner[i]++;
          }

          if ((j >= lfrom) && (j < lto)) {
            m = npartner[j];
            partner[j][m] = tag[i];
            shearpartner[j][m][0] = -shear[0];
            shearpartner[j][m][1] = -shear[1];
            shearpartner[j][m][2] = -shear[2];
            npartner[j]++;
          }
        }
      }
    }

    // set maxtouch = max # of partners of any owned atom
    maxtouch = m = 0;
    for (i = lfrom; i < lto; i++)
      m = MAX(m,npartner[i]);

#if defined(_OPENMP)
#pragma omp critical
#endif
    maxtouch = MAX(m,maxtouch);
  }
}