Added pair style sph/lj/gpu in src/GPU, removed commented codes

lib/gpu/lal_amoeba.cpp

@@ -281,13 +281,7 @@ int AmoebaT::polar_real(const int eflag, const int vflag) {
 
   const int BX=this->block_size();
   const int GX=static_cast<int>(ceil(static_cast<double>(ainum)/(BX/this->_threads_per_atom)));
-  /*
+
-  const int cus = this->device->gpu->cus();
-  while (GX < cus && GX > 1) {
-    BX /= 2;
-    GX=static_cast<int>(ceil(static_cast<double>(ainum)/(BX/this->_threads_per_atom)));
-  }
-  */
   this->time_pair.start();
 
   // Build the short neighbor list if not done yet

lib/gpu/lal_hippo.cpp

@@ -603,13 +603,7 @@ int HippoT::polar_real(const int eflag, const int vflag) {
 
   const int BX=this->block_size();
   const int GX=static_cast<int>(ceil(static_cast<double>(ainum)/(BX/this->_threads_per_atom)));
-  /*
+
-  const int cus = this->device->gpu->cus();
-  while (GX < cus && GX > 1) {
-    BX /= 2;
-    GX=static_cast<int>(ceil(static_cast<double>(ainum)/(BX/this->_threads_per_atom)));
-  }
-  */
   this->time_pair.start();
 
   // Build the short neighbor list if not done yet

lib/gpu/lal_sph_lj.cpp

@@ -107,7 +107,7 @@ int SPHLJT::init(const int ntypes,
     ef_nall=2000;
 
   _max_drhoE_size=static_cast<int>(static_cast<double>(ef_nall)*1.10);
-  drhoE.alloc(_max_drhoE_size,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE);
+  drhoE.alloc(_max_drhoE_size*2,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE);
 
   _allocated=true;
   this->_max_bytes=coeff.row_bytes()+drhoE.row_bytes()+sp_lj.row_bytes();
@@ -134,7 +134,7 @@ double SPHLJT::host_memory_usage() const {
 template <class numtyp, class acctyp>
 void SPHLJT::update_drhoE(void **drhoE_ptr) {
   *drhoE_ptr=drhoE.host.begin();
-  drhoE.update_host(_max_drhoE_size,false);
+  drhoE.update_host(_max_drhoE_size*2,false);
 }
 
 // ---------------------------------------------------------------------------
@@ -148,7 +148,7 @@ int SPHLJT::loop(const int eflag, const int vflag) {
   // Resize drhoE array if necessary
   if (nall > _max_drhoE_size) {
     _max_drhoE_size=static_cast<int>(static_cast<double>(nall)*1.10);
-    drhoE.resize(_max_drhoE_size);
+    drhoE.resize(_max_drhoE_size*2);
   }
 
   // signal that we need to transfer extra data from the host

lib/gpu/lal_sph_lj_ext.cpp

@@ -94,13 +94,15 @@ void sph_lj_gpu_clear() {
 
 int ** sph_lj_gpu_compute_n(const int ago, const int inum_full, const int nall,
                          double **host_x, int *host_type, double *sublo,
-                         double *subhi, tagint *tag, int **nspecial,
+                         double *subhi, int **nspecial,
                          tagint **special, const bool eflag, const bool vflag,
                          const bool eatom, const bool vatom, int &host_start,
                          int **ilist, int **jnum, const double cpu_time, bool &success,
-                         double **host_v, const double dtinvsqrt,
+                         double **host_v, double *boxlo, double *prd) {
-                         const int seed, const int timestep,
+  double dtinvsqrt = 1.0;
-                         double *boxlo, double *prd) {
+  int seed = 0;
+  int timestep = 0;
+  tagint* tag = nullptr;
   return SPHLJMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
                          subhi, tag, nspecial, special, eflag, vflag, eatom,
                          vatom, host_start, ilist, jnum, cpu_time, success,
@@ -108,16 +110,18 @@ int ** sph_lj_gpu_compute_n(const int ago, const int inum_full, const int nall,
 }
 
 void sph_lj_gpu_compute(const int ago, const int inum_full, const int nall,
-                       double **host_x, int *host_type, int *ilist, int *numj,
+                        double **host_x, int *host_type, int *ilist, int *numj,
-                       int **firstneigh, const bool eflag, const bool vflag,
+                        int **firstneigh, const bool eflag, const bool vflag,
-                       const bool eatom, const bool vatom, int &host_start,
+                        const bool eatom, const bool vatom, int &host_start,
-                       const double cpu_time, bool &success, tagint *tag,
+                        const double cpu_time, bool &success,
-                       double **host_v, const double dtinvsqrt,
+                        double **host_v, const int nlocal, double *boxlo, double *prd) {
-                       const int seed, const int timestep,
+  double dtinvsqrt = 1.0;
-                       const int nlocal, double *boxlo, double *prd) {
+  int seed = 0;
+  int timestep = 0;
+  tagint* tag = nullptr;
   SPHLJMF.compute(ago, inum_full, nall, host_x, host_type, ilist, numj,
-                firstneigh, eflag, vflag, eatom, vatom, host_start, cpu_time, success,
+                  firstneigh, eflag, vflag, eatom, vatom, host_start, cpu_time, success,
-                tag, host_v, dtinvsqrt, seed, timestep, nlocal, boxlo, prd);
+                  tag, host_v, dtinvsqrt, seed, timestep, nlocal, boxlo, prd);
 }
 
 void sph_lj_gpu_get_extra_data(double *host_rho, double *host_esph,

src/GPU/pair_sph_lj_gpu.cpp

@@ -0,0 +1,199 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Trung Dac Nguyen (U Chicago)
+------------------------------------------------------------------------- */
+
+#include "pair_sph_lj_gpu.h"
+
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "gpu_extra.h"
+#include "info.h"
+#include "neigh_list.h"
+#include "neighbor.h"
+#include "suffix.h"
+#include "update.h"
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+
+// External functions from cuda library for atom decomposition
+
+int sph_lj_gpu_init(const int ntypes, double **cutsq, double** host_cut,
+                    double **host_viscosity, double *special_lj,
+                    const int inum, const int nall,
+                    const int max_nbors,  const int maxspecial,
+                    const double cell_size, int &gpu_mode, FILE *screen);
+void sph_lj_gpu_clear();
+int **sph_lj_gpu_compute_n(const int ago, const int inum_full, const int nall, double **host_x,
+                        int *host_type, double *sublo, double *subhi, int **nspecial,
+                        tagint **special, const bool eflag, const bool vflag, const bool eatom,
+                        const bool vatom, int &host_start, int **ilist, int **jnum,
+                        const double cpu_time, bool &success, double **host_v,
+                        double *boxlo, double *prd);
+void sph_lj_gpu_compute(const int ago, const int inum_full, const int nall, double **host_x,
+                        int *host_type, int *ilist, int *numj, int **firstneigh, const bool eflag,
+                        const bool vflag, const bool eatom, const bool vatom, int &host_start,
+                        const double cpu_time, bool &success, double **host_v,
+                        const int nlocal, double *boxlo, double *prd);
+void sph_lj_gpu_get_extra_data(double *host_rho, double *host_esph,
+                               double *host_cv, double *host_mass);
+void sph_lj_gpu_update_drhoE(void **drhoE_ptr);
+double sph_lj_gpu_bytes();
+
+/* ---------------------------------------------------------------------- */
+
+PairSPHLJGPU::PairSPHLJGPU(LAMMPS *lmp) : PairSPHLJ(lmp), gpu_mode(GPU_FORCE)
+{
+  drhoE_pinned = nullptr;
+  respa_enable = 0;
+  reinitflag = 0;
+  cpu_time = 0.0;
+  suffix_flag |= Suffix::GPU;
+  GPU_EXTRA::gpu_ready(lmp->modify, lmp->error);
+}
+
+/* ----------------------------------------------------------------------
+   free all arrays
+------------------------------------------------------------------------- */
+
+PairSPHLJGPU::~PairSPHLJGPU()
+{
+  sph_lj_gpu_clear();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSPHLJGPU::compute(int eflag, int vflag)
+{
+  ev_init(eflag, vflag);
+
+  int nall = atom->nlocal + atom->nghost;
+  int inum, host_start;
+
+  bool success = true;
+  int *ilist, *numneigh, **firstneigh;
+
+  double *rho = atom->rho;
+  double *esph = atom->esph;
+  double *cv = atom->cv;
+  double *mass = atom->mass;
+  sph_lj_gpu_get_extra_data(rho, esph, cv, mass);
+
+  if (gpu_mode != GPU_FORCE) {
+    double sublo[3], subhi[3];
+    if (domain->triclinic == 0) {
+      sublo[0] = domain->sublo[0];
+      sublo[1] = domain->sublo[1];
+      sublo[2] = domain->sublo[2];
+      subhi[0] = domain->subhi[0];
+      subhi[1] = domain->subhi[1];
+      subhi[2] = domain->subhi[2];
+    } else {
+      domain->bbox(domain->sublo_lamda, domain->subhi_lamda, sublo, subhi);
+    }
+    inum = atom->nlocal;
+    firstneigh = sph_lj_gpu_compute_n(
+        neighbor->ago, inum, nall, atom->x, atom->type, sublo, subhi, atom->nspecial,
+        atom->special, eflag, vflag, eflag_atom, vflag_atom, host_start, &ilist, &numneigh,
+        cpu_time, success, atom->v, domain->boxlo, domain->prd);
+  } else {
+    inum = list->inum;
+    ilist = list->ilist;
+    numneigh = list->numneigh;
+    firstneigh = list->firstneigh;
+    sph_lj_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type, ilist, numneigh, firstneigh,
+                    eflag, vflag, eflag_atom, vflag_atom, host_start, cpu_time, success,
+                    atom->v, atom->nlocal, domain->boxlo, domain->prd);
+  }
+  if (!success) error->one(FLERR, "Insufficient memory on accelerator");
+
+  // get the drho and dE from device
+
+  double *drho = atom->drho;
+  double *desph = atom->desph;
+  sph_lj_gpu_update_drhoE(&drhoE_pinned);
+
+  int nlocal = atom->nlocal;
+  if (acc_float) {
+    auto drhoE_ptr = (float *)drhoE_pinned;
+    int idx = 0;
+    for (int i = 0; i < nlocal; i++) {
+      drho[i] = drhoE_ptr[idx];
+      desph[i] = drhoE_ptr[idx+1];
+      idx += 2;
+    }
+
+  } else {
+    auto drhoE_ptr = (double *)drhoE_pinned;
+    int idx = 0;
+    for (int i = 0; i < nlocal; i++) {
+      drho[i] = drhoE_ptr[idx];
+      desph[i] = drhoE_ptr[idx+1];
+      idx += 2;
+    }
+  }
+
+  if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)
+    neighbor->build_topology();
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairSPHLJGPU::init_style()
+{
+
+  // Repeat cutsq calculation because done after call to init_style
+  double maxcut = -1.0;
+  double mcut;
+  for (int i = 1; i <= atom->ntypes; i++) {
+    for (int j = i; j <= atom->ntypes; j++) {
+      if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) {
+        mcut = init_one(i, j);
+        mcut *= mcut;
+        if (mcut > maxcut) maxcut = mcut;
+        cutsq[i][j] = cutsq[j][i] = mcut;
+      } else
+        cutsq[i][j] = cutsq[j][i] = 0.0;
+    }
+  }
+  double cell_size = sqrt(maxcut) + neighbor->skin;
+
+  int maxspecial = 0;
+  if (atom->molecular != Atom::ATOMIC) maxspecial = atom->maxspecial;
+  int mnf = 5e-2 * neighbor->oneatom;
+  int success =
+      sph_lj_gpu_init(atom->ntypes + 1, cutsq, cut, viscosity,
+                      force->special_lj, atom->nlocal, atom->nlocal + atom->nghost,
+                      mnf, maxspecial, cell_size, gpu_mode, screen);
+  GPU_EXTRA::check_flag(success, error, world);
+
+  acc_float = Info::has_accelerator_feature("GPU", "precision", "single");
+
+  if (gpu_mode == GPU_FORCE) neighbor->add_request(this, NeighConst::REQ_FULL);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairSPHLJGPU::memory_usage()
+{
+  double bytes = Pair::memory_usage();
+  return bytes + sph_lj_gpu_bytes();
+}

src/GPU/pair_sph_lj_gpu.h

@@ -0,0 +1,48 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   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.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(sph/lj/gpu,PairSPHLJGPU);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_SPH_LJ_GPU_H
+#define LMP_PAIR_SPH_LJ_GPU_H
+
+#include "pair_sph_lj.h"
+
+namespace LAMMPS_NS {
+
+class PairSPHLJGPU : public PairSPHLJ {
+ public:
+  PairSPHLJGPU(LAMMPS *lmp);
+  ~PairSPHLJGPU() override;
+  void cpu_compute(int, int, int, int, int *, int *, int **);
+  void compute(int, int) override;
+  void init_style() override;
+  double memory_usage() override;
+
+  enum { GPU_FORCE, GPU_NEIGH, GPU_HYB_NEIGH };
+
+  void *drhoE_pinned;
+  bool acc_float;
+
+ private:
+  int gpu_mode;
+  double cpu_time;
+};
+
+}    // namespace LAMMPS_NS
+#endif
+#endif