ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Adding sph/taitwater in lib/gpu

Browse Source
This commit is contained in:
Trung Nguyen
2023-12-07 17:22:45 -06:00
parent 6fe16c7606
commit cad4c25750
8 changed files with 835 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/gpu/lal_edpd.cpp
View File

@ -104,7 +104,7 @@ int EDPDT::init(const int ntypes,
this->atom->type_pack4(ntypes,lj_types,coeff2,host_write,host_power,host_kappa, this->atom->type_pack4(ntypes,lj_types,coeff2,host_write,host_power,host_kappa,
host_powerT,host_cutT); host_powerT,host_cutT);
UCL_H_Vec<int> dview_mass(ntypes, *(this->ucl_device), UCL_WRITE_ONLY); UCL_H_Vec<numtyp> dview_mass(ntypes, *(this->ucl_device), UCL_WRITE_ONLY);
for (int i = 0; i < ntypes; i++) for (int i = 0; i < ntypes; i++)
dview_mass[i] = host_mass[i]; dview_mass[i] = host_mass[i];
mass.alloc(ntypes,*(this->ucl_device), UCL_READ_ONLY); mass.alloc(ntypes,*(this->ucl_device), UCL_READ_ONLY);

4
lib/gpu/lal_sph_lj.cpp
View File

@ -163,8 +163,8 @@ int SPHLJT::loop(const int eflag, const int vflag) {
int idx = n+i*nstride; int idx = n+i*nstride;
numtyp4 v; numtyp4 v;
v.x = rho[i]; v.x = rho[i];
v.x = esph[i]; v.y = esph[i];
v.y = cv[i]; v.z = cv[i];
v.w = mass[i]; v.w = mass[i];
pextra[idx] = v; pextra[idx] = v;
} }

4
lib/gpu/lal_sph_lj.cu
View File

@ -1,9 +1,9 @@
// ************************************************************************** // **************************************************************************
// edpd.cu // sph_lj.cu
// ------------------- // -------------------
// Trung Dac Nguyen (U Chicago) // Trung Dac Nguyen (U Chicago)
// //
// Device code for acceleration of the edpd pair style // Device code for acceleration of the sph/lj pair style
// //
// __________________________________________________________________________ // __________________________________________________________________________
// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL) // This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)

221
lib/gpu/lal_sph_taitwater.cpp Normal file
View File

@ -0,0 +1,221 @@
/***************************************************************************
sph_taitwater.cpp
-------------------
Trung Dac Nguyen (U Chicago)
Class for acceleration of the sph_taitwater pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : December 2023
email : ndactrung@gmail.com
***************************************************************************/
#if defined(USE_OPENCL)
#include "sph_taitwater_cl.h"
#elif defined(USE_CUDART)
const char *sph_taitwater=0;
#else
#include "sph_taitwater_cubin.h"
#endif
#include "lal_sph_taitwater.h"
#include <cassert>
namespace LAMMPS_AL {
#define SPHTaitwaterT SPHTaitwater<numtyp, acctyp>
extern Device<PRECISION,ACC_PRECISION> device;
template <class numtyp, class acctyp>
SPHTaitwaterT::SPHTaitwater() : BaseDPD<numtyp,acctyp>(), _allocated(false) {
_max_drhoE_size = 0;
}
template <class numtyp, class acctyp>
SPHTaitwaterT::~SPHTaitwater() {
clear();
}
template <class numtyp, class acctyp>
int SPHTaitwaterT::bytes_per_atom(const int max_nbors) const {
return this->bytes_per_atom_atomic(max_nbors);
}
template <class numtyp, class acctyp>
int SPHTaitwaterT::init(const int ntypes,
double **host_cutsq, double **host_cut,
double **host_viscosity, double* host_rho0,
double* host_soundspeed, double *host_special_lj,
const int nlocal, const int nall,
const int max_nbors, const int maxspecial,
const double cell_size,
const double gpu_split, FILE *_screen) {
const int max_shared_types=this->device->max_shared_types();
int onetype=0;
#ifdef USE_OPENCL
if (maxspecial==0)
for (int i=1; i<ntypes; i++)
for (int j=i; j<ntypes; j++)
if (host_cutsq[i][j]>0) {
if (onetype>0)
onetype=-1;
else if (onetype==0)
onetype=i*max_shared_types+j;
}
if (onetype<0) onetype=0;
#endif
int success;
int extra_fields = 4; // round up to accomodate quadruples of numtyp values
// rho, mass
success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,
gpu_split,_screen,sph_taitwater,"k_sph_taitwater",onetype,extra_fields);
if (success!=0)
return success;
// If atom type constants fit in shared memory use fast kernel
int lj_types=ntypes;
shared_types=false;
if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
lj_types=max_shared_types;
shared_types=true;
}
_lj_types=taitwater_types;
// Allocate a host write buffer for data initialization
UCL_H_Vec<numtyp> host_write(lj_types*lj_types*32,*(this->ucl_device),
UCL_WRITE_ONLY);
for (int i=0; i<lj_types*lj_types; i++)
host_write[i]=0.0;
coeff.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,coeff,host_write,host_viscosity,
host_cut, host_cutsq);
UCL_H_Vec<numtyp2> dview_rho0ss(ntypes, *(this->ucl_device), UCL_WRITE_ONLY);
for (int i = 0; i < ntypes; i++) {
dview_rho0ss[i].x = host_rho0[i];
dview_rho0ss[i].y = host_soundspeed[i];
}
rho0sspeed.alloc(ntypes,*(this->ucl_device), UCL_READ_ONLY);
ucl_copy(rho0sspeed,dview_rho0ss,false);
UCL_H_Vec<double> dview;
sp_lj.alloc(4,*(this->ucl_device),UCL_READ_ONLY);
dview.view(host_special_lj,4,*(this->ucl_device));
ucl_copy(sp_lj,dview,false);
// allocate per-atom array Q
int ef_nall=nall;
if (ef_nall==0)
ef_nall=2000;
_max_drhoE_size=static_cast<int>(static_cast<double>(ef_nall)*1.10);
drhoE.alloc(_max_drhoE_size*2,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE);
_allocated=true;
this->_max_bytes=coeff.row_bytes()+corho0sspeedeff.row_bytes()+drhoE.row_bytes()+sp_lj.row_bytes();
return 0;
}
template <class numtyp, class acctyp>
void SPHTaitwaterT::clear() {
if (!_allocated)
return;
_allocated=false;
coeff.clear();
dview_rho0ss.clear();
drhoE.clear();
sp_lj.clear();
this->clear_atomic();
}
template <class numtyp, class acctyp>
double SPHTaitwaterT::host_memory_usage() const {
return this->host_memory_usage_atomic()+sizeof(SPHTaitwater<numtyp,acctyp>);
}
template <class numtyp, class acctyp>
void SPHTaitwaterT::update_drhoE(void **drhoE_ptr) {
*drhoE_ptr=drhoE.host.begin();
drhoE.update_host(_max_drhoE_size*2,false);
}
// ---------------------------------------------------------------------------
// Calculate energies, forces, and torques
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
int SPHTaitwaterT::loop(const int eflag, const int vflag) {
int nall = this->atom->nall();
// 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*2);
}
// signal that we need to transfer extra data from the host
this->atom->extra_data_unavail();
numtyp4 *pextra=reinterpret_cast<numtyp4*>(&(this->atom->extra[0]));
int n = 0;
int nstride = 1;
for (int i = 0; i < nall; i++) {
int idx = n+i*nstride;
numtyp4 v;
v.x = rho[i];
v.y = mass[i];
v.z = 0;
v.w = 0;
pextra[idx] = v;
}
this->atom->add_extra_data();
// Compute the block size and grid size to keep all cores busy
const int BX=this->block_size();
int GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
(BX/this->_threads_per_atom)));
int ainum=this->ans->inum();
int nbor_pitch=this->nbor->nbor_pitch();
this->time_pair.start();
if (shared_types) {
this->k_pair_sel->set_size(GX,BX);
this->k_pair_sel->run(&this->atom->x, &this->atom->extra, &coeff, &rho0sspeed, &sp_lj,
&this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &drhoE, &eflag, &vflag,
&ainum, &nbor_pitch, &this->atom->v, &this->_threads_per_atom);
} else {
this->k_pair.set_size(GX,BX);
this->k_pair.run(&this->atom->x, &this->atom->extra, &coeff, &rho0sspeed,
&_lj_types, &sp_lj, &this->nbor->dev_nbor, &this->_nbor_data->begin(),
&this->ans->force, &this->ans->engv, &drhoE, &eflag, &vflag,
&ainum, &nbor_pitch, &this->atom->v, &this->_threads_per_atom);
}
this->time_pair.stop();
return GX;
}
// ---------------------------------------------------------------------------
// Get the extra data pointers from host
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void SPHTaitwaterT::get_extra_data(double *host_rho, double* host_mass) {
rho = host_rho;
mass = host_mass;
}
template class SPHTaitwater<PRECISION,ACC_PRECISION>;
}

374
lib/gpu/lal_sph_taitwater.cu Normal file
View File

@ -0,0 +1,374 @@
// **************************************************************************
// sph_taitwater.cu
// -------------------
// Trung Dac Nguyen (U Chicago)
//
// Device code for acceleration of the sph/taitwater pair style
//
// __________________________________________________________________________
// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
// __________________________________________________________________________
//
// begin : September 2023
// email : ndactrung@gmail.com
// ***************************************************************************
#if defined(NV_KERNEL) || defined(USE_HIP)
#include "lal_aux_fun1.h"
#ifndef _DOUBLE_DOUBLE
_texture( pos_tex,float4);
_texture( vel_tex,float4);
#else
_texture_2d( pos_tex,int4);
_texture_2d( vel_tex,int4);
#endif
#else
#define pos_tex x_
#define vel_tex v_
#endif
#if (SHUFFLE_AVAIL == 0)
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE) \
if (t_per_atom>1) { \
simdsync(); \
simd_reduce_add2(t_per_atom, red_acc, offset, tid, \
drhoEacc.x, drhoEacc.y); \
} \
if (offset==0 && ii<inum) { \
drhoE[ii]=drhoEacc; \
}
#else
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE) \
if (t_per_atom>1) { \
for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \
drhoEacc.x += shfl_down(drhoEacc.x, s, t_per_atom); \
drhoEacc.y += shfl_down(drhoEacc.y, s, t_per_atom); \
} \
} \
if (offset==0 && ii<inum) { \
drhoE[ii].x=drhoEacc; \
}
#endif
__kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict extra,
const __global numtyp4 *restrict coeff,
const __global numtyp2 *restrict rho0sspeed,
const int lj_types,
const __global numtyp *restrict sp_lj,
const __global int * dev_nbor,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int t_per_atom) {
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
int n_stride;
local_allocate_store_pair();
acctyp3 f;
f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
acctyp energy, virial[6];
if (EVFLAG) {
energy=(acctyp)0;
for (int i=0; i<6; i++) virial[i]=(acctyp)0;
}
acctyp2 drhoEacc;
drhoEacc.x = drhoEacc.x = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int itype=ix.w;
numtyp4 iv; fetch4(iv,i,vel_tex); //v_[i];
const numtyp4 extrai = extra[i];
numtyp rhoi = extrai.x;
numtyp massi= extrai.y;
// compute pressure of atom i with Tait EOS
numtyp tmp = rho[i] / rho0[itype];
numtyp fi = tmp * tmp * tmp;
fi = B[itype] * (fi * fi * tmp - 1.0);
fi /= (rhoi * rhoi);
for ( ; nbor<nbor_end; nbor+=n_stride) {
ucl_prefetch(dev_packed+nbor+n_stride);
int j=dev_packed[nbor];
factor_lj = sp_lj[sbmask(j)];
j &= NEIGHMASK;
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
int jtype=jx.w;
numtyp4 jv; fetch4(jv,j,vel_tex); //v_[j];
// Compute r12
numtyp delx = ix.x-jx.x;
numtyp dely = ix.y-jx.y;
numtyp delz = ix.z-jx.z;
numtyp rsq = delx*delx+dely*dely+delz*delz;
int mtype=itype*lj_types+jtype;
if (rsq<coeff[mtype].z) { // cutsq[itype][jtype]
const numtyp coeffx=coeff[mtype].x; // viscosity[itype][jtype]
const numtyp coeffy=coeff[mtype].y; // cut[itype][jtype]
const numtyp4 extraj = extra[j];
numtyp rhoj = extraj.x;
numtyp esphj = extraj.y;
numtyp cvj = extraj.z;
numtyp massj= extraj.w;
numtyp h = coeffy; // cut[itype][jtype]
ih = ucl_recip(h); // (numtyp)1.0 / h;
numtyp ihsq = ih * ih;
numtyp ihcub = ihsq * ih;
numtyp wfd = h - ucl_sqrt(rsq);
// domain->dimension == 3 Lucy Kernel, 3d
wfd = (numtyp)-25.066903536973515383 * wfd * wfd * ihsq * ihsq * ihsq * ih;
// Lucy Kernel, 2d
//wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq;
// compute pressure of atom j with Tait EOS
numtyp tmp = rho[j] / rho0[jtype];
numtyp fj = tmp * tmp * tmp;
fj = B[jtype] * (fj * fj * tmp - 1.0);
fj /= (rhoj * rhoj);
// dot product of velocity delta and distance vector
numtyp delvx = iv.x - jv.x;
numtyp delvy = iv.y - jv.y;
numtyp delvz = iv.z - jv.z;
numtyp delVdotDelR = delx*delvx + dely*delvy + delz*delvz;
// artificial viscosity (Monaghan 1992)
numtyp fvisc = (numtyp)0;
if (delVdotDelR < (numyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numyp)0.01 * h * h);
//fvisc = -coeffx * (ci + cj) * mu / (rhoi + rhoj); // viscosity[itype][jtype]
fvisc = -coeffx * (soundspeed[itype]
+ soundspeed[jtype]) * mu / (rhoi + rhoj);
}
// total pair force & thermal energy increment
numtyp force = -massi * massj * (fi + fj + fvisc) * wfd;
numtyp deltaE = (numtyp)-0.5 * force * delVdotDelR;
f.x+=delx*force;
f.y+=dely*force;
f.z+=delz*force;
// and change in density, drho[i]
drhoEacc.x += massj * delVdotDelR * wfd;
// change in thermal energy, desph[i]
drhoEacc.y += deltaE;
if (EVFLAG && eflag) {
numtyp e = (numtyp)0;
energy+=e;
}
if (EVFLAG && vflag) {
virial[0] += delx*delx*force;
virial[1] += dely*dely*force;
virial[2] += delz*delz*force;
virial[3] += delx*dely*force;
virial[4] += delx*delz*force;
virial[5] += dely*delz*force;
}
}
} // for nbor
} // if ii
store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag,
ans,engv);
store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,drhoE);
}
__kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict extra,
const __global numtyp4 *restrict coeff_in,
const __global numtyp2 *restrict rho0sspeed_in,
const __global numtyp *restrict sp_lj_in,
const __global int * dev_nbor,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int t_per_atom) {
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
#ifndef ONETYPE
__local numtyp4 coeff[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
__local numtyp sp_lj[4];
if (tid<4) {
sp_lj[tid]=sp_lj_in[tid];
}
if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
coeff[tid]=coeff_in[tid];
}
__syncthreads();
#else
const numtyp coeffx=coeff_in[ONETYPE].x; // viscosity[itype][jtype]
const numtyp coeffy=coeff_in[ONETYPE].y; // cut[itype][jtype]
const numtyp cutsq_p=coeff_in[ONETYPE].z; // cutsq[itype][jtype]
#endif
int n_stride;
local_allocate_store_pair();
acctyp3 f;
f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
acctyp energy, virial[6];
if (EVFLAG) {
energy=(acctyp)0;
for (int i=0; i<6; i++) virial[i]=(acctyp)0;
}
acctyp2 drhoEacc;
drhoEacc.x = drhoEacc.x = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int iw=ix.w;
numtyp mass_itype = mass[iw];
#ifndef ONETYPE
int itype=fast_mul((int)MAX_SHARED_TYPES,iw);
#endif
numtyp4 iv; fetch4(iv,i,vel_tex); //v_[i];
int itag=iv.w;
const numtyp4 extrai = extra[i];
numtyp rhoi = extrai.x;
numtyp esphi = extrai.y;
numtyp cvi = extrai.z;
numtyp massi= extrai.w;
// compute pressure of particle i with LJ EOS
numtyp fci[2];
LJEOS2(rhoi, esphi, cvi, fci);
numtyp fi = fci[0];
numtyp ci = fci[1];
fi /= (rhoi * rhoi);
for ( ; nbor<nbor_end; nbor+=n_stride) {
ucl_prefetch(dev_packed+nbor+n_stride);
int j=dev_packed[nbor];
#ifndef ONETYPE
factor_dpd = sp_lj[sbmask(j)];
j &= NEIGHMASK;
#endif
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
#ifndef ONETYPE
int mtype=itype+jx.w;
const numtyp cutsq_p=cutsq[mtype];
#endif
numtyp4 jv; fetch4(jv,j,vel_tex); //v_[j];
int jtag=jv.w;
// Compute r12
numtyp delx = ix.x-jx.x;
numtyp dely = ix.y-jx.y;
numtyp delz = ix.z-jx.z;
numtyp rsq = delx*delx+dely*dely+delz*delz;
if (rsq<cutsq_p) {
#ifndef ONETYPE
const numtyp coeffx=coeff[mtype].x; // viscosity[itype][jtype]
const numtyp coeffy=coeff[mtype].y; // cut[itype][jtype]
#endif
const numtyp4 extraj = extra[j];
numtyp rhoj = extraj.x;
numtyp esphj = extraj.y;
numtyp cvj = extraj.z;
numtyp massj= extraj.w;
numtyp h = coeffy; // cut[itype][jtype]
ih = ih = ucl_recip(h); // (numtyp)1.0 / h;
numtyp ihsq = ih * ih;
numtyp ihcub = ihsq * ih;
numtyp wfd = h - ucl_sqrt(rsq);
// domain->dimension == 3 Lucy Kernel, 3d
wfd = (numtyp)-25.066903536973515383 * wfd * wfd * ihsq * ihsq * ihsq * ih;
// Lucy Kernel, 2d
//wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq;
// function call to LJ EOS
numtyp fcj[2];
LJEOS2(rhoj, esphj, cvj, fcj);
numtyp fj = fcj[0];
numtyp cj = fcj[1];
fj /= (rhoj * rhoj);
// apply long-range correction to model a LJ fluid with cutoff
// this implies that the modelled LJ fluid has cutoff == SPH cutoff
numtyp lrc = (numtyp)-11.1701 * (ihcub * ihcub * ihcub - (numtyp)1.5 * ihcub);
fi += lrc;
fj += lrc;
// dot product of velocity delta and distance vector
numtyp delvx = iv.x - jv.x;
numtyp delvy = iv.y - jv.y;
numtyp delvz = iv.z - jv.z;
numtyp delVdotDelR = delx*delvx + dely*delvy + delz*delvz;
// artificial viscosity (Monaghan 1992)
numtyp fvisc = (numtyp)0;
if (delVdotDelR < (numyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numyp)0.01 * h * h);
fvisc = -coeffx * (ci + cj) * mu / (rhoi + rhoj); // viscosity[itype][jtype]
}
// total pair force & thermal energy increment
numtyp force = -massi * massj * (fi + fj + fvisc) * wfd;
numtyp deltaE = (numtyp)-0.5 * force * delVdotDelR;
f.x+=delx*force;
f.y+=dely*force;
f.z+=delz*force;
if (EVFLAG && eflag) {
numtyp e = (numtyp)0;
energy+=e;
}
if (EVFLAG && vflag) {
virial[0] += delx*delx*force;
virial[1] += dely*dely*force;
virial[2] += delz*delz*force;
virial[3] += delx*dely*force;
virial[4] += delx*delz*force;
virial[5] += dely*delz*force;
}
}
} // for nbor
} // if ii
store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag, ans,engv);
store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,drhoE);
}

93
lib/gpu/lal_sph_taitwater.h Normal file
View File

@ -0,0 +1,93 @@
/***************************************************************************
sph_lj.h
-------------------
Trung Dac Nguyen (U Chicago)
Class for acceleration of the sph lj pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : December 2023
email : ndactrung@gmail.com
***************************************************************************/
#ifndef LAL_SPH_TAITWATER_H
#define LAL_SPH_TaitLAL_SPH_TAITWATER_Hwater_H
#include "lal_base_dpd.h"
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class SPHTaitwater : public BaseDPD<numtyp, acctyp> {
public:
SPHTaitwater();
~SPHTaitwater();
/// Clear any previous data and set up for a new LAMMPS run
/** \param max_nbors initial number of rows in the neighbor matrix
* \param cell_size cutoff + skin
* \param gpu_split fraction of particles handled by device
*
* Returns:
* - 0 if successful
* - -1 if fix gpu not found
* - -3 if there is an out of memory error
* - -4 if the GPU library was not compiled for GPU
* - -5 Double precision is not supported on card **/
int init(const int ntypes, double **host_cutsq,
double** host_cut, double **host_viscosity,
double* host_rho0, double* host_soundspeed,
double *host_special_lj, const int nlocal,
const int nall, const int max_nbors,
const int maxspecial, const double cell_size,
const double gpu_split, FILE *screen);
/// Clear all host and device data
/** \note This is called at the beginning of the init() routine **/
void clear();
/// Returns memory usage on device per atom
int bytes_per_atom(const int max_nbors) const;
/// Total host memory used by library for pair style
double host_memory_usage() const;
void get_extra_data(double *host_rho, double* host_mass);
/// copy drho and desph from device to host
void update_drhoE(void **drhoE_ptr);
// --------------------------- TYPE DATA --------------------------
/// coeff.x = viscosity, coeff.y = cut, coeff.z = cutsq
UCL_D_Vec<numtyp2> coeff;
UCL_D_Vec<numtyp2> rho0sspeed;
/// Special LJ values
UCL_D_Vec<numtyp> sp_lj;
/// If atom type constants fit in shared memory, use fast kernels
bool shared_types;
/// Number of atom types
int _lj_types;
/// Per-atom arrays
UCL_Vector<acctyp,acctyp> drhoE;
int _max_drhoE_size;
/// pointer to host data
double *rho, *mass;
private:
bool _allocated;
int loop(const int eflag, const int vflag);
};
}
#endif

140
lib/gpu/lal_sph_taitwater_ext.cpp Normal file
View File

@ -0,0 +1,140 @@
/***************************************************************************
sph_taitwater_ext.cpp
-------------------
Trung Dac Nguyen (U Chicago)
Functions for LAMMPS access to sph lj acceleration routines.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : December 2023
email : ndactrung@gmail.com
***************************************************************************/
#include <iostream>
#include <cassert>
#include <cmath>
#include "lal_sph_taitwater.h"
using namespace std;
using namespace LAMMPS_AL;
static SPHTaitwater<PRECISION,ACC_PRECISION> SPHTaitwaterMF;
// ---------------------------------------------------------------------------
// Allocate memory on host and device and copy constants to device
// ---------------------------------------------------------------------------
int sph_taitwater_gpu_init(const int ntypes, double **cutsq, double** host_cut,
double **host_viscosity, double* host_rho0,
double* host_soundspeed, 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) {
SPHTaitwaterMF.clear();
gpu_mode=SPHTaitwaterMF.device->gpu_mode();
double gpu_split=SPHTaitwaterMF.device->particle_split();
int first_gpu=SPHTaitwaterMF.device->first_device();
int last_gpu=SPHTaitwaterMF.device->last_device();
int world_me=SPHTaitwaterMF.device->world_me();
int gpu_rank=SPHTaitwaterMF.device->gpu_rank();
int procs_per_gpu=SPHTaitwaterMF.device->procs_per_gpu();
SPHTaitwaterMF.device->init_message(screen,"sph_taitwater",first_gpu,last_gpu);
bool message=false;
if (SPHTaitwaterMF.device->replica_me()==0 && screen)
message=true;
if (message) {
fprintf(screen,"Initializing Device and compiling on process 0...");
fflush(screen);
}
int init_ok=0;
if (world_me==0)
init_ok=SPHTaitwaterMF.init(ntypes, cutsq, host_cut, host_viscosity,
host_rho0, host_soundspeed, special_lj,
inum, nall, max_nbors, maxspecial, cell_size,
gpu_split, screen);
SPHTaitwaterMF.device->world_barrier();
if (message)
fprintf(screen,"Done.\n");
for (int i=0; i<procs_per_gpu; i++) {
if (message) {
if (last_gpu-first_gpu==0)
fprintf(screen,"Initializing Device %d on core %d...",first_gpu,i);
else
fprintf(screen,"Initializing Devices %d-%d on core %d...",first_gpu,
last_gpu,i);
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=SPHTaitwaterMF.init(ntypes, cutsq, host_cut, host_viscosity,
host_rho0, host_soundspeed, special_lj,
inum, nall, max_nbors, maxspecial, cell_size,
gpu_split, screen);
SPHTaitwaterMF.device->serialize_init();
if (message)
fprintf(screen,"Done.\n");
}
if (message)
fprintf(screen,"\n");
if (init_ok==0)
SPHTaitwaterMF.estimate_gpu_overhead();
return init_ok;
}
void sph_taitwater_gpu_clear() {
SPHTaitwaterMF.clear();
}
int ** sph_taitwater_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) {
double dtinvsqrt = 1.0;
int seed = 0;
int timestep = 0;
tagint* tag = nullptr;
return SPHTaitwaterMF.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,
host_v, dtinvsqrt, seed, timestep, boxlo, prd);
}
void sph_taitwater_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) {
double dtinvsqrt = 1.0;
int seed = 0;
int timestep = 0;
tagint* tag = nullptr;
SPHTaitwaterMF.compute(ago, inum_full, nall, host_x, host_type, ilist, numj,
firstneigh, eflag, vflag, eatom, vatom, host_start, cpu_time, success,
tag, host_v, dtinvsqrt, seed, timestep, nlocal, boxlo, prd);
}
void sph_taitwater_gpu_get_extra_data(double *host_rho, double *host_mass) {
SPHTaitwaterMF.get_extra_data(host_rho, host_mass);
}
void sph_taitwater_gpu_update_drhoE(void **drhoE_ptr) {
SPHTaitwaterMF.update_drhoE(drhoE_ptr);
}
double sph_taitwater_gpu_bytes() {
return SPHTaitwaterMF.host_memory_usage();
}

4
src/GPU/pair_sph_lj_gpu.cpp
View File

@ -118,8 +118,8 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
numneigh = list->numneigh; numneigh = list->numneigh;
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
sph_lj_gpu_compute(neighbor->ago, inum, nall, atom->x, atom->type, ilist, numneigh, 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, eflag, vflag, eflag_atom, vflag_atom, host_start, cpu_time, success,
atom->v, atom->nlocal, domain->boxlo, domain->prd); atom->v, atom->nlocal, domain->boxlo, domain->prd);
} }
if (!success) error->one(FLERR, "Insufficient memory on accelerator"); if (!success) error->one(FLERR, "Insufficient memory on accelerator");