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

Added BaseSPH for sph pair styles, fixed bugs with per-type arrays

Browse Source
This commit is contained in:
Trung Nguyen
2023-12-08 14:09:39 -06:00
parent 0f78afc66f
commit 46a670889d
12 changed files with 750 additions and 165 deletions
Download Patch File Download Diff File Expand all files Collapse all files

362
lib/gpu/lal_base_sph.cpp Normal file
View File

@ -0,0 +1,362 @@
/***************************************************************************
base_sph.cpp
-------------------
Trung Dac Nguyen (ORNL)
Base class for SPH pair styles needing per-particle data for position,
velocity, and type.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : December 2023
email : ndactrung@gmail.com
***************************************************************************/
#include "lal_base_sph.h"
namespace LAMMPS_AL {
#define BaseSPHT BaseSPH<numtyp, acctyp>
extern Device<PRECISION,ACC_PRECISION> global_device;
template <class numtyp, class acctyp>
BaseSPHT::BaseSPH() : _compiled(false), _max_bytes(0) {
device=&global_device;
ans=new Answer<numtyp,acctyp>();
nbor=new Neighbor();
pair_program=nullptr;
ucl_device=nullptr;
#if defined(LAL_OCL_EV_JIT)
pair_program_noev=nullptr;
#endif
}
template <class numtyp, class acctyp>
BaseSPHT::~BaseSPH() {
delete ans;
delete nbor;
k_pair_fast.clear();
k_pair.clear();
if (pair_program) delete pair_program;
#if defined(LAL_OCL_EV_JIT)
k_pair_noev.clear();
if (pair_program_noev) delete pair_program_noev;
#endif
}
template <class numtyp, class acctyp>
int BaseSPHT::bytes_per_atom_atomic(const int max_nbors) const {
return device->atom.bytes_per_atom()+ans->bytes_per_atom()+
nbor->bytes_per_atom(max_nbors);
}
template <class numtyp, class acctyp>
int BaseSPHT::init_atomic(const int nlocal, const int nall,
const int max_nbors, const int maxspecial,
const double cell_size, const double gpu_split,
FILE *_screen, const void *pair_program,
const char *k_name, const int onetype,
const int extra_fields) {
screen=_screen;
int gpu_nbor=0;
if (device->gpu_mode()==Device<numtyp,acctyp>::GPU_NEIGH)
gpu_nbor=1;
else if (device->gpu_mode()==Device<numtyp,acctyp>::GPU_HYB_NEIGH)
gpu_nbor=2;
int _gpu_host=0;
int host_nlocal=hd_balancer.first_host_count(nlocal,gpu_split,gpu_nbor);
if (host_nlocal>0)
_gpu_host=1;
_threads_per_atom=device->threads_per_atom();
bool charge = false;
bool rot = false;
bool vel = true;
_extra_fields = extra_fields;
int success=device->init(*ans,charge,rot,nlocal,nall,maxspecial,vel,_extra_fields/4);
if (success!=0)
return success;
if (ucl_device!=device->gpu) _compiled=false;
ucl_device=device->gpu;
atom=&device->atom;
_block_size=device->pair_block_size();
compile_kernels(*ucl_device,pair_program,k_name,onetype);
if (_threads_per_atom>1 && gpu_nbor==0) {
nbor->packing(true);
_nbor_data=&(nbor->dev_packed);
} else
_nbor_data=&(nbor->dev_nbor);
success = device->init_nbor(nbor,nlocal,host_nlocal,nall,maxspecial,_gpu_host,
max_nbors,cell_size,false,_threads_per_atom);
if (success!=0)
return success;
// Initialize host-device load balancer
hd_balancer.init(device,gpu_nbor,gpu_split);
// Initialize timers for the selected GPU
time_pair.init(*ucl_device);
time_pair.zero();
pos_tex.bind_float(atom->x,4);
vel_tex.bind_float(atom->v,4);
_max_an_bytes=ans->gpu_bytes()+nbor->gpu_bytes();
return success;
}
template <class numtyp, class acctyp>
void BaseSPHT::estimate_gpu_overhead() {
device->estimate_gpu_overhead(1,_gpu_overhead,_driver_overhead);
}
template <class numtyp, class acctyp>
void BaseSPHT::clear_atomic() {
// Output any timing information
acc_timers();
double avg_split=hd_balancer.all_avg_split();
_gpu_overhead*=hd_balancer.timestep();
_driver_overhead*=hd_balancer.timestep();
device->output_times(time_pair,*ans,*nbor,avg_split,_max_bytes+_max_an_bytes,
_gpu_overhead,_driver_overhead,_threads_per_atom,screen);
time_pair.clear();
hd_balancer.clear();
nbor->clear();
ans->clear();
}
// ---------------------------------------------------------------------------
// Copy neighbor list from host
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
int * BaseSPHT::reset_nbors(const int nall, const int inum, int *ilist,
int *numj, int **firstneigh, bool &success) {
success=true;
int mn=nbor->max_nbor_loop(inum,numj,ilist);
resize_atom(inum,nall,success);
resize_local(inum,mn,success);
if (!success)
return nullptr;
nbor->get_host(inum,ilist,numj,firstneigh,block_size());
double bytes=ans->gpu_bytes()+nbor->gpu_bytes();
if (bytes>_max_an_bytes)
_max_an_bytes=bytes;
return ilist;
}
// ---------------------------------------------------------------------------
// Build neighbor list on device
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
inline void BaseSPHT::build_nbor_list(const int inum, const int host_inum,
const int nall, double **host_x,
int *host_type, double *sublo,
double *subhi, tagint *tag,
int **nspecial, tagint **special,
bool &success) {
success=true;
resize_atom(inum,nall,success);
resize_local(inum,host_inum,nbor->max_nbors(),success);
if (!success)
return;
atom->cast_copy_x(host_x,host_type);
int mn;
nbor->build_nbor_list(host_x, inum, host_inum, nall, *atom, sublo, subhi,
tag, nspecial, special, success, mn, ans->error_flag);
double bytes=ans->gpu_bytes()+nbor->gpu_bytes();
if (bytes>_max_an_bytes)
_max_an_bytes=bytes;
}
// ---------------------------------------------------------------------------
// Copy nbor list from host if necessary and then calculate forces, virials,..
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void BaseSPHT::compute(const int f_ago, const int inum_full, const int nall,
double **host_x, int *host_type, int *ilist, int *numj,
int **firstneigh, const bool eflag_in, const bool vflag_in,
const bool eatom, const bool vatom, int &host_start,
const double cpu_time, bool &success, tagint *tag,
double **host_v, const int nlocal) {
acc_timers();
int eflag, vflag;
if (eatom) eflag=2;
else if (eflag_in) eflag=1;
else eflag=0;
if (vatom) vflag=2;
else if (vflag_in) vflag=1;
else vflag=0;
#ifdef LAL_NO_BLOCK_REDUCE
if (eflag) eflag=2;
if (vflag) vflag=2;
#endif
set_kernel(eflag,vflag);
if (inum_full==0) {
host_start=0;
// Make sure textures are correct if realloc by a different hybrid style
resize_atom(0,nall,success);
zero_timers();
return;
}
int ago=hd_balancer.ago_first(f_ago);
int inum=hd_balancer.balance(ago,inum_full,cpu_time);
ans->inum(inum);
host_start=inum;
if (ago==0) {
reset_nbors(nall, inum, ilist, numj, firstneigh, success);
if (!success)
return;
}
atom->cast_x_data(host_x,host_type);
atom->cast_v_data(host_v,tag);
hd_balancer.start_timer();
atom->add_x_data(host_x,host_type);
atom->add_v_data(host_v,tag);
const int red_blocks=loop(eflag,vflag);
ans->copy_answers(eflag_in,vflag_in,eatom,vatom,ilist,red_blocks);
device->add_ans_object(ans);
hd_balancer.stop_timer();
}
// ---------------------------------------------------------------------------
// Reneighbor on GPU if necessary and then compute forces, virials, energies
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
int** BaseSPHT::compute(const int ago, const int inum_full, const int nall,
double **host_x, int *host_type, double *sublo,
double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag_in, const bool vflag_in,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **jnum, const double cpu_time, bool &success,
double **host_v) {
acc_timers();
int eflag, vflag;
if (eatom) eflag=2;
else if (eflag_in) eflag=1;
else eflag=0;
if (vatom) vflag=2;
else if (vflag_in) vflag=1;
else vflag=0;
#ifdef LAL_NO_BLOCK_REDUCE
if (eflag) eflag=2;
if (vflag) vflag=2;
#endif
set_kernel(eflag,vflag);
if (inum_full==0) {
host_start=0;
// Make sure textures are correct if realloc by a different hybrid style
resize_atom(0,nall,success);
zero_timers();
return nullptr;
}
hd_balancer.balance(cpu_time);
int inum=hd_balancer.get_gpu_count(ago,inum_full);
ans->inum(inum);
host_start=inum;
// Build neighbor list on GPU if necessary
if (ago==0) {
build_nbor_list(inum, inum_full-inum, nall, host_x, host_type,
sublo, subhi, tag, nspecial, special, success);
if (!success)
return nullptr;
atom->cast_v_data(host_v,tag);
hd_balancer.start_timer();
} else {
atom->cast_x_data(host_x,host_type);
atom->cast_v_data(host_v,tag);
hd_balancer.start_timer();
atom->add_x_data(host_x,host_type);
}
atom->add_v_data(host_v,tag);
*ilist=nbor->host_ilist.begin();
*jnum=nbor->host_acc.begin();
const int red_blocks=loop(eflag,vflag);
ans->copy_answers(eflag_in,vflag_in,eatom,vatom,red_blocks);
device->add_ans_object(ans);
hd_balancer.stop_timer();
return nbor->host_jlist.begin()-host_start;
}
template <class numtyp, class acctyp>
double BaseSPHT::host_memory_usage_atomic() const {
return device->atom.host_memory_usage()+nbor->host_memory_usage()+
4*sizeof(numtyp)+sizeof(BaseSPH<numtyp,acctyp>);
}
template <class numtyp, class acctyp>
void BaseSPHT::compile_kernels(UCL_Device &dev, const void *pair_str,
const char *kname, const int onetype) {
if (_compiled && _onetype==onetype)
return;
_onetype=onetype;
std::string s_fast=std::string(kname)+"_fast";
if (pair_program) delete pair_program;
pair_program=new UCL_Program(dev);
std::string oclstring = device->compile_string()+" -DEVFLAG=1";
if (_onetype) oclstring+=" -DONETYPE="+device->toa(_onetype);
pair_program->load_string(pair_str,oclstring.c_str(),nullptr,screen);
k_pair_fast.set_function(*pair_program,s_fast.c_str());
k_pair.set_function(*pair_program,kname);
pos_tex.get_texture(*pair_program,"pos_tex");
vel_tex.get_texture(*pair_program,"vel_tex");
#if defined(LAL_OCL_EV_JIT)
oclstring = device->compile_string()+" -DEVFLAG=0";
if (_onetype) oclstring+=" -DONETYPE="+device->toa(_onetype);
if (pair_program_noev) delete pair_program_noev;
pair_program_noev=new UCL_Program(dev);
pair_program_noev->load_string(pair_str,oclstring.c_str(),nullptr,screen);
k_pair_noev.set_function(*pair_program_noev,s_fast.c_str());
#else
k_pair_sel = &k_pair_fast;
#endif
_compiled=true;
#if defined(USE_OPENCL) && (defined(CL_VERSION_2_1) || defined(CL_VERSION_3_0))
if (dev.has_subgroup_support()) {
size_t mx_subgroup_sz = k_pair_fast.max_subgroup_size(_block_size);
#if defined(LAL_OCL_EV_JIT)
mx_subgroup_sz = std::min(mx_subgroup_sz, k_pair_noev.max_subgroup_size(_block_size));
#endif
if (_threads_per_atom > (int)mx_subgroup_sz) _threads_per_atom = mx_subgroup_sz;
device->set_simd_size(mx_subgroup_sz);
}
#endif
}
template class BaseSPH<PRECISION,ACC_PRECISION>;
}

209
lib/gpu/lal_base_sph.h Normal file
View File

@ -0,0 +1,209 @@
/***************************************************************************
base_sph.h
-------------------
Trung Dac Nguyen (U Chicago)
Base class for SPH pair styles needing per-particle data for position,
velocity, and type.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
__________________________________________________________________________
begin : December 2023
email : ndactrung@gmail.com
***************************************************************************/
#ifndef LAL_BASE_SPH_H
#define LAL_BASE_DPD_H
#include "lal_device.h"
#include "lal_balance.h"
#include "mpi.h"
#ifdef USE_OPENCL
#include "geryon/ocl_texture.h"
#elif defined(USE_HIP)
#include "geryon/hip_texture.h"
#else
#include "geryon/nvd_texture.h"
#endif
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class BaseSPH {
public:
BaseSPH();
virtual ~BaseSPH();
/// 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
* \param k_name name for the kernel for force calculation
*
* 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_atomic(const int nlocal, const int nall, const int max_nbors,
const int maxspecial, const double cell_size,
const double gpu_split, FILE *screen,
const void *pair_program, const char *k_name,
const int onetype=0, const int extra_fields=0);
/// Estimate the overhead for GPU context changes and CPU driver
void estimate_gpu_overhead();
/// Check if there is enough storage for atom arrays and realloc if not
/** \param success set to false if insufficient memory **/
inline void resize_atom(const int inum, const int nall, bool &success) {
if (atom->resize(nall, success)) {
pos_tex.bind_float(atom->x,4);
vel_tex.bind_float(atom->v,4);
}
ans->resize(inum,success);
}
/// Check if there is enough storage for neighbors and realloc if not
/** \param nlocal number of particles whose nbors must be stored on device
* \param host_inum number of particles whose nbors need to copied to host
* \param current maximum number of neighbors
* \note olist_size=total number of local particles **/
inline void resize_local(const int inum, const int max_nbors, bool &success) {
nbor->resize(inum,max_nbors,success);
}
/// Check if there is enough storage for neighbors and realloc if not
/** \param nlocal number of particles whose nbors must be stored on device
* \param host_inum number of particles whose nbors need to copied to host
* \param current maximum number of neighbors
* \note host_inum is 0 if the host is performing neighboring
* \note nlocal+host_inum=total number local particles
* \note olist_size=0 **/
inline void resize_local(const int inum, const int host_inum,
const int max_nbors, bool &success) {
nbor->resize(inum,host_inum,max_nbors,success);
}
/// Clear all host and device data
/** \note This is called at the beginning of the init() routine **/
void clear_atomic();
/// Returns memory usage on device per atom
int bytes_per_atom_atomic(const int max_nbors) const;
/// Total host memory used by library for pair style
double host_memory_usage_atomic() const;
/// Accumulate timers
inline void acc_timers() {
if (device->time_device()) {
nbor->acc_timers(screen);
time_pair.add_to_total();
atom->acc_timers();
ans->acc_timers();
}
}
/// Zero timers
inline void zero_timers() {
time_pair.zero();
atom->zero_timers();
ans->zero_timers();
}
/// Copy neighbor list from host
int * reset_nbors(const int nall, const int inum, int *ilist, int *numj,
int **firstneigh, bool &success);
/// Build neighbor list on device
void build_nbor_list(const int inum, const int host_inum,
const int nall, double **host_x, int *host_type,
double *sublo, double *subhi, tagint *tag, int **nspecial,
tagint **special, bool &success);
/// Pair loop with host neighboring
void compute(const int f_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, tagint *tag,
double **v, const int nlocal);
/// Pair loop with device neighboring
int** compute(const int ago, const int inum_full, const int nall,
double **host_x, int *host_type, double *sublo,
double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **numj, const double cpu_time, bool &success,
double **v);
// -------------------------- DEVICE DATA -------------------------
/// Device Properties and Atom and Neighbor storage
Device<numtyp,acctyp> *device;
/// Geryon device
UCL_Device *ucl_device;
/// Device Timers
UCL_Timer time_pair;
/// Host device load balancer
Balance<numtyp,acctyp> hd_balancer;
/// LAMMPS pointer for screen output
FILE *screen;
// --------------------------- ATOM DATA --------------------------
/// Atom Data
Atom<numtyp,acctyp> *atom;
// ------------------------ FORCE/ENERGY DATA -----------------------
Answer<numtyp,acctyp> *ans;
// --------------------------- NBOR DATA ----------------------------
/// Neighbor data
Neighbor *nbor;
// ------------------------- DEVICE KERNELS -------------------------
UCL_Program *pair_program, *pair_program_noev;
UCL_Kernel k_pair_fast, k_pair, k_pair_noev, *k_pair_sel;
inline int block_size() { return _block_size; }
inline void set_kernel(const int eflag, const int vflag) {
#if defined(LAL_OCL_EV_JIT)
if (eflag || vflag) k_pair_sel = &k_pair_fast;
else k_pair_sel = &k_pair_noev;
#endif
}
// --------------------------- TEXTURES -----------------------------
UCL_Texture pos_tex;
UCL_Texture vel_tex;
// ------------------------- COMMON VARS ----------------------------
protected:
bool _compiled;
int _block_size, _threads_per_atom, _onetype, _extra_fields;
double _max_bytes, _max_an_bytes;
double _gpu_overhead, _driver_overhead;
UCL_D_Vec<int> *_nbor_data;
void compile_kernels(UCL_Device &dev, const void *pair_string,
const char *k, const int onetype);
virtual int loop(const int eflag, const int vflag) = 0;
};
}
#endif

25
lib/gpu/lal_sph_lj.cpp
View File

@ -29,7 +29,7 @@ namespace LAMMPS_AL {
extern Device<PRECISION,ACC_PRECISION> device;
template <class numtyp, class acctyp>
SPHLJT::SPHLJ() : BaseDPD<numtyp,acctyp>(), _allocated(false) {
SPHLJT::SPHLJ() : BaseSPH<numtyp,acctyp>(), _allocated(false) {
_max_drhoE_size = 0;
}
@ -46,8 +46,8 @@ int SPHLJT::bytes_per_atom(const int max_nbors) const {
template <class numtyp, class acctyp>
int SPHLJT::init(const int ntypes,
double **host_cutsq, double **host_cut,
double **host_viscosity, const int dimension,
double *host_special_lj,
double **host_viscosity, double* host_mass,
const int dimension, double *host_special_lj,
const int nlocal, const int nall,
const int max_nbors, const int maxspecial,
const double cell_size,
@ -70,7 +70,7 @@ int SPHLJT::init(const int ntypes,
int success;
int extra_fields = 4; // round up to accomodate quadruples of numtyp values
// rho, cv, mass
// rho, cv
success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,
gpu_split,_screen,sph_lj,"k_sph_lj",onetype,extra_fields);
if (success!=0)
@ -96,6 +96,12 @@ int SPHLJT::init(const int ntypes,
this->atom->type_pack4(ntypes,lj_types,coeff,host_write,host_viscosity,
host_cut, host_cutsq);
UCL_H_Vec<numtyp> dview_mass(ntypes, *(this->ucl_device), UCL_WRITE_ONLY);
for (int i = 0; i < ntypes; i++)
dview_mass[i] = host_mass[i];
mass.alloc(ntypes,*(this->ucl_device), UCL_READ_ONLY);
ucl_copy(mass,dview_mass,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));
@ -124,6 +130,7 @@ void SPHLJT::clear() {
_allocated=false;
coeff.clear();
mass.clear();
drhoE.clear();
sp_lj.clear();
this->clear_atomic();
@ -168,7 +175,7 @@ int SPHLJT::loop(const int eflag, const int vflag) {
v.x = rho[i];
v.y = esph[i];
v.z = cv[i];
v.w = mass[i];
v.w = 0;
pextra[idx] = v;
}
this->atom->add_extra_data();
@ -184,13 +191,13 @@ int SPHLJT::loop(const int eflag, const int vflag) {
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, &sp_lj,
this->k_pair_sel->run(&this->atom->x, &this->atom->extra, &coeff, &mass, &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, &_dimension, &this->_threads_per_atom);
} else {
this->k_pair.set_size(GX,BX);
this->k_pair.run(&this->atom->x, &this->atom->extra, &coeff,
this->k_pair.run(&this->atom->x, &this->atom->extra, &coeff, &mass,
&_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, &_dimension, &this->_threads_per_atom);
@ -205,12 +212,10 @@ int SPHLJT::loop(const int eflag, const int vflag) {
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void SPHLJT::get_extra_data(double *host_rho, double *host_esph,
double *host_cv, double* host_mass) {
void SPHLJT::get_extra_data(double *host_rho, double *host_esph, double *host_cv) {
rho = host_rho;
esph = host_esph;
cv = host_cv;
mass = host_mass;
}
template class SPHLJ<PRECISION,ACC_PRECISION>;

27
lib/gpu/lal_sph_lj.cu
View File

@ -47,7 +47,7 @@ _texture_2d( vel_tex,int4);
} \
} \
if (offset==0 && ii<inum) { \
drhoE[ii].x=drhoEacc; \
drhoE[ii]=drhoEacc; \
}
#endif
@ -98,13 +98,14 @@ void LJEOS2(const numtyp rho, const numtyp e, const numtyp cv, numtyp pc[2])
__kernel void k_sph_lj(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict extra,
const __global numtyp4 *restrict coeff,
const __global numtyp *restrict mass,
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,
__global acctyp2 *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
@ -132,13 +133,13 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int itype=ix.w;
numtyp mass_itype = mass[itype];
numtyp4 iv; fetch4(iv,i,vel_tex); //v_[i];
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];
@ -151,7 +152,6 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
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];
@ -166,7 +166,7 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
int mtype=itype*lj_types+jtype;
if (rsq<coeff[mtype].z) { // cutsq[itype][jtype]
numtyp mass_jtype = mass[jtype];
const numtyp coeffx=coeff[mtype].x; // viscosity[itype][jtype]
const numtyp coeffy=coeff[mtype].y; // cut[itype][jtype]
@ -174,7 +174,6 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
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;
@ -187,7 +186,7 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
wfd = (numtyp)-25.066903536973515383 * wfd * wfd * ihsq * ihsq * ihsq * ih;
} else {
// Lucy Kernel, 2d
wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq;
wfd = (numtyp)-19.098593171027440292 * wfd * wfd * ihsq * ihsq * ihsq;
}
// function call to LJ EOS
@ -217,7 +216,7 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
}
// total pair force & thermal energy increment
numtyp force = -massi * massj * (fi + fj + fvisc) * wfd;
numtyp force = -mass_itype * mass_type * (fi + fj + fvisc) * wfd;
numtyp deltaE = (numtyp)-0.5 * force * delVdotDelR;
f.x+=delx*force;
@ -253,12 +252,13 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
__kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict extra,
const __global numtyp4 *restrict coeff_in,
const __global numtyp *restrict mass,
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,
__global acctyp2 *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
@ -327,7 +327,6 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
int j=dev_packed[nbor];
#ifndef ONETYPE
factor_dpd = sp_lj[sbmask(j)];
j &= NEIGHMASK;
#endif
@ -367,7 +366,7 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
wfd = (numtyp)-25.066903536973515383 * wfd * wfd * ihsq * ihsq * ihsq * ih;
} else {
// Lucy Kernel, 2d
wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq;
wfd = (numtyp)-19.098593171027440292 * wfd * wfd * ihsq * ihsq * ihsq;
}
// function call to LJ EOS
@ -404,6 +403,12 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
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;

16
lib/gpu/lal_sph_lj.h
View File

@ -16,12 +16,12 @@
#ifndef LAL_SPH_LJ_H
#define LAL_SPH_LJ_H
#include "lal_base_dpd.h"
#include "lal_base_sph.h"
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class SPHLJ : public BaseDPD<numtyp, acctyp> {
class SPHLJ : public BaseSPH<numtyp, acctyp> {
public:
SPHLJ();
~SPHLJ();
@ -38,7 +38,8 @@ class SPHLJ : public BaseDPD<numtyp, acctyp> {
* - -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, const int dimension,
double** host_cut, double **host_viscosity, double *host_mass,
const int dimension,
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);
@ -54,7 +55,7 @@ class SPHLJ : public BaseDPD<numtyp, acctyp> {
double host_memory_usage() const;
void get_extra_data(double *host_rho, double *host_esph,
double *host_cv, double* host_mass);
double *host_cv);
/// copy drho and desph from device to host
void update_drhoE(void **drhoE_ptr);
@ -62,7 +63,10 @@ class SPHLJ : public BaseDPD<numtyp, acctyp> {
// --------------------------- TYPE DATA --------------------------
/// coeff.x = viscosity, coeff.y = cut, coeff.z = cutsq
UCL_D_Vec<numtyp2> coeff;
UCL_D_Vec<numtyp4> coeff;
/// per-type coeffs
UCL_D_Vec<numtyp> mass;
/// Special LJ values
UCL_D_Vec<numtyp> sp_lj;
@ -80,7 +84,7 @@ class SPHLJ : public BaseDPD<numtyp, acctyp> {
int _dimension;
/// pointer to host data
double *rho, *esph, *cv, *mass;
double *rho, *esph, *cv;
private:
bool _allocated;

45
lib/gpu/lal_sph_lj_ext.cpp
View File

@ -28,7 +28,7 @@ static SPHLJ<PRECISION,ACC_PRECISION> SPHLJMF;
// Allocate memory on host and device and copy constants to device
// ---------------------------------------------------------------------------
int sph_lj_gpu_init(const int ntypes, double **cutsq, double** host_cut,
double **host_viscosity, const int dimension,
double **host_viscosity, double* host_mass, const int dimension,
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) {
@ -54,8 +54,8 @@ int sph_lj_gpu_init(const int ntypes, double **cutsq, double** host_cut,
int init_ok=0;
if (world_me==0)
init_ok=SPHLJMF.init(ntypes, cutsq, host_cut, host_viscosity, dimension,
special_lj, inum, nall, max_nbors, maxspecial,
init_ok=SPHLJMF.init(ntypes, cutsq, host_cut, host_viscosity, host_mass,
dimension, special_lj, inum, nall, max_nbors, maxspecial,
cell_size, gpu_split, screen);
SPHLJMF.device->world_barrier();
@ -72,8 +72,8 @@ int sph_lj_gpu_init(const int ntypes, double **cutsq, double** host_cut,
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=SPHLJMF.init(ntypes, cutsq, host_cut, host_viscosity, dimension,
special_lj, inum, nall, max_nbors, maxspecial,
init_ok=SPHLJMF.init(ntypes, cutsq, host_cut, host_viscosity, host_mass,
dimension, special_lj, inum, nall, max_nbors, maxspecial,
cell_size, gpu_split, screen);
SPHLJMF.device->serialize_init();
@ -93,38 +93,31 @@ 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,
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;
double **host_x, int *host_type, double *sublo,
double *subhi, tagint *host_tag, 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) {
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,
host_v, dtinvsqrt, seed, timestep, boxlo, prd);
subhi, host_tag, nspecial, special, eflag, vflag,
eatom, vatom, host_start, ilist, jnum, cpu_time, success,
host_v);
}
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, tagint *tag,
double **host_v, const int nlocal, double *boxlo, double *prd) {
double dtinvsqrt = 1.0;
int seed = 0;
int timestep = 0;
const double cpu_time, bool &success, tagint *host_tag,
double **host_v, const int nlocal) {
SPHLJMF.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);
host_tag, host_v, nlocal);
}
void sph_lj_gpu_get_extra_data(double *host_rho, double *host_esph,
double *host_cv, double *host_mass) {
SPHLJMF.get_extra_data(host_rho, host_esph, host_cv, host_mass);
void sph_lj_gpu_get_extra_data(double *host_rho, double *host_esph, double *host_cv) {
SPHLJMF.get_extra_data(host_rho, host_esph, host_cv);
}
void sph_lj_gpu_update_drhoE(void **drhoE_ptr) {

25
lib/gpu/lal_sph_taitwater.cpp
View File

@ -3,7 +3,7 @@
-------------------
Trung Dac Nguyen (U Chicago)
Class for acceleration of the sph_taitwater pair style.
Class for acceleration of the sph/taitwater pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
@ -29,7 +29,7 @@ namespace LAMMPS_AL {
extern Device<PRECISION,ACC_PRECISION> device;
template <class numtyp, class acctyp>
SPHTaitwaterT::SPHTaitwater() : BaseDPD<numtyp,acctyp>(), _allocated(false) {
SPHTaitwaterT::SPHTaitwater() : BaseSPH<numtyp,acctyp>(), _allocated(false) {
_max_drhoE_size = 0;
}
@ -46,8 +46,8 @@ int SPHTaitwaterT::bytes_per_atom(const int max_nbors) const {
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_B, const int dimension,
double* host_mass, double* host_rho0,
double* host_soundspeed, double* host_B, const int dimension,
double *host_special_lj, const int nlocal, const int nall,
const int max_nbors, const int maxspecial,
const double cell_size,
@ -70,7 +70,7 @@ int SPHTaitwaterT::init(const int ntypes, double **host_cutsq,
int success;
int extra_fields = 4; // round up to accomodate quadruples of numtyp values
// rho, mass
// rho
success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,
gpu_split,_screen,sph_taitwater,"k_sph_taitwater",
onetype,extra_fields);
@ -99,10 +99,10 @@ int SPHTaitwaterT::init(const int ntypes, double **host_cutsq,
UCL_H_Vec<numtyp4> dview_coeff2(ntypes, *(this->ucl_device), UCL_WRITE_ONLY);
for (int i = 0; i < ntypes; i++) {
dview_coeff2[i].x = host_rho0[i];
dview_coeff2[i].y = host_soundspeed[i];
dview_coeff2[i].z = host_B[i];
dview_coeff2[i].w = 0;
dview_coeff2[i].x = host_mass[i];
dview_coeff2[i].y = host_rho0[i];
dview_coeff2[i].z = host_soundspeed[i];
dview_coeff2[i].w = host_B[i];
}
coeff2.alloc(ntypes,*(this->ucl_device), UCL_READ_ONLY);
ucl_copy(coeff2,dview_coeff2,false);
@ -122,7 +122,7 @@ int SPHTaitwaterT::init(const int ntypes, double **host_cutsq,
drhoE.alloc(_max_drhoE_size*2,*(this->ucl_device),UCL_READ_WRITE,UCL_READ_WRITE);
_dimension = dimension;
_allocated=true;
this->_max_bytes=coeff.row_bytes()+coeff2.row_bytes()+drhoE.row_bytes()+sp_lj.row_bytes();
return 0;
@ -178,7 +178,7 @@ int SPHTaitwaterT::loop(const int eflag, const int vflag) {
int idx = n+i*nstride;
numtyp4 v;
v.x = rho[i];
v.y = mass[i];
v.y = 0;
v.z = 0;
v.w = 0;
pextra[idx] = v;
@ -217,9 +217,8 @@ int SPHTaitwaterT::loop(const int eflag, const int vflag) {
// ---------------------------------------------------------------------------
template <class numtyp, class acctyp>
void SPHTaitwaterT::get_extra_data(double *host_rho, double* host_mass) {
void SPHTaitwaterT::get_extra_data(double *host_rho) {
rho = host_rho;
mass = host_mass;
}
template class SPHTaitwater<PRECISION,ACC_PRECISION>;

74
lib/gpu/lal_sph_taitwater.cu
View File

@ -47,7 +47,7 @@ _texture_2d( vel_tex,int4);
} \
} \
if (offset==0 && ii<inum) { \
drhoE[ii].x=drhoEacc; \
drhoE[ii]=drhoEacc; \
}
#endif
@ -61,7 +61,7 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp *restrict drhoE,
__global acctyp2 *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
@ -89,14 +89,14 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int itype=ix.w;
numtyp rho0_itype = coeff2[itype].x;
numtyp soundspeed_itype = coeff2[itype].y;
numtyp B_itype = coeff2[itype].z;
numtyp mass_itype = coeff2[itype].x;
numtyp rho0_itype = coeff2[itype].y;
numtyp soundspeed_itype = coeff2[itype].z;
numtyp B_itype = coeff2[itype].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 = rhoi / rho0_itype;
@ -108,7 +108,6 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
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];
@ -126,16 +125,16 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
const numtyp coeffx=coeff[mtype].x; // viscosity[itype][jtype]
const numtyp coeffy=coeff[mtype].y; // cut[itype][jtype]
numtyp rho0_jtype = coeff2[jtype].x;
numtyp soundspeed_jtype = coeff2[jtype].y;
numtyp B_jtype = coeff2[jtype].z;
numtyp mass_jtype = coeff2[jtype].x;
numtyp rho0_jtype = coeff2[jtype].y;
numtyp soundspeed_jtype = coeff2[jtype].z;
numtyp B_jtype = coeff2[jtype].w;
const numtyp4 extraj = extra[j];
numtyp rhoj = extraj.x;
numtyp massj= extraj.y;
numtyp h = coeffy; // cut[itype][jtype]
ih = ucl_recip(h); // (numtyp)1.0 / h;
numtyp ih = ucl_recip(h); // (numtyp)1.0 / h;
numtyp ihsq = ih * ih;
numtyp ihcub = ihsq * ih;
@ -145,7 +144,7 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
wfd = (numtyp)-25.066903536973515383 * wfd * wfd * ihsq * ihsq * ihsq * ih;
} else {
// Lucy Kernel, 2d
wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq;
wfd = (numtyp)-19.098593171027440292 * wfd * wfd * ihsq * ihsq * ihsq;
}
// compute pressure of atom j with Tait EOS
@ -163,14 +162,14 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
// artificial viscosity (Monaghan 1992)
numtyp fvisc = (numtyp)0;
if (delVdotDelR < (numyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numyp)0.01 * h * h);
if (delVdotDelR < (numtyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numtyp)0.01 * h * h);
fvisc = -coeffx * (soundspeed_itype
+ soundspeed_jtype) * mu / (rhoi + rhoj);
}
// total pair force & thermal energy increment
numtyp force = -massi * massj * (fi + fj + fvisc) * wfd;
numtyp force = -mass_itype * mass_jtype * (fi + fj + fvisc) * wfd;
numtyp deltaE = (numtyp)-0.5 * force * delVdotDelR;
f.x+=delx*force;
@ -178,7 +177,7 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
f.z+=delz*force;
// and change in density, drho[i]
drhoEacc.x += massj * delVdotDelR * wfd;
drhoEacc.x += mass_jtype* delVdotDelR * wfd;
// change in thermal energy, desph[i]
drhoEacc.y += deltaE;
@ -207,13 +206,12 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
const __global numtyp4 *restrict extra,
const __global numtyp4 *restrict coeff_in,
const __global numtyp4 *restrict coeff2_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,
__global acctyp2 *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
@ -223,10 +221,14 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
#ifndef ONETYPE
__local numtyp4 coeff[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
__local numtyp4 coeff2[MAX_SHARED_TYPES];
__local numtyp sp_lj[4];
if (tid<4) {
sp_lj[tid]=sp_lj_in[tid];
}
if (tid<MAX_SHARED_TYPES) {
coeff2[tid] = coeff2_in[tid];
}
if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
coeff[tid]=coeff_in[tid];
}
@ -257,9 +259,10 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int iw=ix.w;
numtyp rho0_itype = coeff2[iw].x;
numtyp soundspeed_itype = coeff2[iw].y;
numtyp B_itype = coeff2[iw].z;
numtyp mass_itype = coeff2[iw].x;
numtyp rho0_itype = coeff2[iw].y;
numtyp soundspeed_itype = coeff2[iw].z;
numtyp B_itype = coeff2[iw].w;
#ifndef ONETYPE
int itype=fast_mul((int)MAX_SHARED_TYPES,iw);
#endif
@ -281,14 +284,14 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
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];
int jtype=jx.w;
#ifndef ONETYPE
int mtype=itype+jx.w;
const numtyp cutsq_p=cutsq[mtype];
const numtyp cutsq_p=coeff[mtype].z;
#endif
numtyp4 jv; fetch4(jv,j,vel_tex); //v_[j];
int jtag=jv.w;
@ -305,16 +308,17 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
const numtyp coeffy=coeff[mtype].y; // cut[itype][jtype]
#endif
numtyp rho0_jtype = coeff2[jtype].x;
numtyp soundspeed_jtype = coeff2[jtype].y;
numtyp B_jtype = coeff2[jtype].z;
numtyp mass_jtype = coeff2[jtype].x;
numtyp rho0_jtype = coeff2[jtype].y;
numtyp soundspeed_jtype = coeff2[jtype].z;
numtyp B_jtype = coeff2[jtype].w;
const numtyp4 extraj = extra[j];
numtyp rhoj = extraj.x;
numtyp massj= extraj.y;
numtyp h = coeffy; // cut[itype][jtype]
ih = ih = ucl_recip(h); // (numtyp)1.0 / h;
numtyp ih = ucl_recip(h); // (numtyp)1.0 / h;
numtyp ihsq = ih * ih;
numtyp ihcub = ihsq * ih;
@ -324,7 +328,7 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
wfd = (numtyp)-25.066903536973515383 * wfd * wfd * ihsq * ihsq * ihsq * ih;
} else {
// Lucy Kernel, 2d
wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq;
wfd = (numtyp)-19.098593171027440292 * wfd * wfd * ihsq * ihsq * ihsq;
}
// compute pressure of atom j with Tait EOS
@ -341,20 +345,26 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
// artificial viscosity (Monaghan 1992)
numtyp fvisc = (numtyp)0;
if (delVdotDelR < (numyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numyp)0.01 * h * h);
if (delVdotDelR < (numtyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numtyp)0.01 * h * h);
fvisc = -coeffx * (soundspeed_itype
+ soundspeed_jtype) * mu / (rhoi + rhoj);
}
// total pair force & thermal energy increment
numtyp force = -massi * massj * (fi + fj + fvisc) * wfd;
numtyp force = -mass_itype * mass_jtype * (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
drhoEacc.x += mass_jtype * delVdotDelR * wfd;
// change in thermal energy
drhoEacc.y += deltaE;
if (EVFLAG && eflag) {
numtyp e = (numtyp)0;
energy+=e;

20
lib/gpu/lal_sph_taitwater.h
View File

@ -1,9 +1,9 @@
/***************************************************************************
sph_lj.h
sph_taitwater.h
-------------------
Trung Dac Nguyen (U Chicago)
Class for acceleration of the sph lj pair style.
Class for acceleration of the sph/taitwater pair style.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
@ -14,14 +14,14 @@
***************************************************************************/
#ifndef LAL_SPH_TAITWATER_H
#define LAL_SPH_TaitLAL_SPH_TAITWATER_Hwater_H
#define LAL_SPH_TAITWATER_H
#include "lal_base_dpd.h"
#include "lal_base_sph.h"
namespace LAMMPS_AL {
template <class numtyp, class acctyp>
class SPHTaitwater : public BaseDPD<numtyp, acctyp> {
class SPHTaitwater : public BaseSPH<numtyp, acctyp> {
public:
SPHTaitwater();
~SPHTaitwater();
@ -38,9 +38,9 @@ class SPHTaitwater : public BaseDPD<numtyp, acctyp> {
* - -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_B, const int dimension, double *host_special_lj,
double** host_cut, double **host_viscosity, double *host_mass,
double* host_rho0, double* host_soundspeed, double* host_B,
const int dimension, 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);
@ -55,7 +55,7 @@ class SPHTaitwater : public BaseDPD<numtyp, acctyp> {
/// Total host memory used by library for pair style
double host_memory_usage() const;
void get_extra_data(double *host_rho, double* host_mass);
void get_extra_data(double *host_rho);
/// copy drho and desph from device to host
void update_drhoE(void **drhoE_ptr);
@ -84,7 +84,7 @@ class SPHTaitwater : public BaseDPD<numtyp, acctyp> {
int _dimension;
/// pointer to host data
double *rho, *mass;
double *rho;
private:
bool _allocated;

34
lib/gpu/lal_sph_taitwater_ext.cpp
View File

@ -3,7 +3,7 @@
-------------------
Trung Dac Nguyen (U Chicago)
Functions for LAMMPS access to sph lj acceleration routines.
Functions for LAMMPS access to sph taitwater acceleration routines.
__________________________________________________________________________
This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
@ -28,8 +28,8 @@ 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* host_B,
double **host_viscosity, double* host_mass,
double* host_rho0, double* host_soundspeed, double* host_B,
const int dimension, double *special_lj,
const int inum, const int nall,
const int max_nbors, const int maxspecial,
@ -56,7 +56,7 @@ int sph_taitwater_gpu_init(const int ntypes, double **cutsq, double** host_cut,
int init_ok=0;
if (world_me==0)
init_ok=SPHTaitwaterMF.init(ntypes, cutsq, host_cut, host_viscosity,
init_ok=SPHTaitwaterMF.init(ntypes, cutsq, host_cut, host_viscosity, host_mass,
host_rho0, host_soundspeed, host_B, dimension,
special_lj, inum, nall, max_nbors, maxspecial,
cell_size, gpu_split, screen);
@ -75,7 +75,7 @@ int sph_taitwater_gpu_init(const int ntypes, double **cutsq, double** host_cut,
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=SPHTaitwaterMF.init(ntypes, cutsq, host_cut, host_viscosity,
init_ok=SPHTaitwaterMF.init(ntypes, cutsq, host_cut, host_viscosity, host_mass,
host_rho0, host_soundspeed, host_B, dimension,
special_lj, inum, nall, max_nbors, maxspecial,
cell_size, gpu_split, screen);
@ -98,36 +98,30 @@ void sph_taitwater_gpu_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, tagint *tag, int **nspecial,
double *subhi, tagint *host_tag, 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;
double **host_v) {
return SPHTaitwaterMF.compute(ago, inum_full, nall, host_x, host_type, sublo,
subhi, tag, nspecial, special, eflag, vflag, eatom,
subhi, host_tag, nspecial, special, eflag, vflag, eatom,
vatom, host_start, ilist, jnum, cpu_time, success,
host_v, dtinvsqrt, seed, timestep, boxlo, prd);
host_v);
}
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, tagint *tag,
double **host_v, const int nlocal, double *boxlo, double *prd) {
double dtinvsqrt = 1.0;
int seed = 0;
int timestep = 0;
const double cpu_time, bool &success, tagint *host_tag,
double **host_v, const int nlocal) {
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);
host_tag, host_v, nlocal);
}
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_get_extra_data(double *host_rho) {
SPHTaitwaterMF.get_extra_data(host_rho);
}
void sph_taitwater_gpu_update_drhoE(void **drhoE_ptr) {

41
src/GPU/pair_sph_lj_gpu.cpp
View File

@ -35,24 +35,27 @@ 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, const int dimension, double *special_lj,
double **host_viscosity, double* host_mass,
const int dimension, 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, tagint *tag, 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, tagint *tag, double **host_v,
const int nlocal, double *boxlo, double *prd);
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 *host_tag, 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);
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, tagint *host_tag,
double **host_v, const int nlocal);
void sph_lj_gpu_get_extra_data(double *host_rho, double *host_esph,
double *host_cv, double *host_mass);
double *host_cv);
void sph_lj_gpu_update_drhoE(void **drhoE_ptr);
double sph_lj_gpu_bytes();
@ -92,8 +95,7 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
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);
sph_lj_gpu_get_extra_data(rho, esph, cv);
if (gpu_mode != GPU_FORCE) {
double sublo[3], subhi[3];
@ -112,7 +114,7 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
neighbor->ago, inum, nall, atom->x, atom->type,
sublo, subhi, atom->tag, atom->nspecial, atom->special, eflag, vflag,
eflag_atom, vflag_atom, host_start, &ilist, &numneigh,
cpu_time, success, atom->v, domain->boxlo, domain->prd);
cpu_time, success, atom->v);
} else {
inum = list->inum;
ilist = list->ilist;
@ -121,7 +123,7 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
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->tag, atom->v, atom->nlocal, domain->boxlo, domain->prd);
atom->tag, atom->v, atom->nlocal);
}
if (!success) error->one(FLERR, "Insufficient memory on accelerator");
@ -182,8 +184,9 @@ void PairSPHLJGPU::init_style()
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, domain->dimension,
force->special_lj, atom->nlocal, atom->nlocal + atom->nghost,
sph_lj_gpu_init(atom->ntypes + 1, cutsq, cut, viscosity, atom->mass,
domain->dimension, force->special_lj, atom->nlocal,
atom->nlocal + atom->nghost,
mnf, maxspecial, cell_size, gpu_mode, screen);
GPU_EXTRA::check_flag(success, error, world);

37
src/GPU/pair_sph_taitwater_gpu.cpp
View File

@ -35,24 +35,26 @@ using namespace LAMMPS_NS;
// External functions from cuda library for atom decomposition
int sph_taitwater_gpu_init(const int ntypes, double **cutsq, double** host_cut,
double **host_viscosity, double* host_rho0,
double **host_viscosity, double* host_mass, double* host_rho0,
double* host_soundspeed, double* host_B, const int dimension,
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_taitwater_gpu_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, tagint *tag, 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_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, tagint *tag, double **host_v,
const int nlocal, double *boxlo, double *prd);
void sph_taitwater_gpu_get_extra_data(double *host_rho, double *host_mass);
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, tagint *tag, 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);
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, tagint *tag,
double **host_v, const int nlocal);
void sph_taitwater_gpu_get_extra_data(double *host_rho);
void sph_taitwater_gpu_update_drhoE(void **drhoE_ptr);
double sph_taitwater_gpu_bytes();
@ -90,8 +92,7 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
int *ilist, *numneigh, **firstneigh;
double *rho = atom->rho;
double *mass = atom->mass;
sph_taitwater_gpu_get_extra_data(rho, mass);
sph_taitwater_gpu_get_extra_data(rho);
if (gpu_mode != GPU_FORCE) {
double sublo[3], subhi[3];
@ -109,7 +110,7 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
firstneigh = sph_taitwater_gpu_compute_n(
neighbor->ago, inum, nall, atom->x, atom->type, sublo, subhi, atom->tag, atom->nspecial,
atom->special, eflag, vflag, eflag_atom, vflag_atom, host_start, &ilist, &numneigh,
cpu_time, success, atom->v, domain->boxlo, domain->prd);
cpu_time, success, atom->v);
} else {
inum = list->inum;
ilist = list->ilist;
@ -117,7 +118,7 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
firstneigh = list->firstneigh;
sph_taitwater_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->tag, atom->v, atom->nlocal, domain->boxlo, domain->prd);
atom->tag, atom->v, atom->nlocal);
}
if (!success) error->one(FLERR, "Insufficient memory on accelerator");
@ -178,7 +179,7 @@ void PairSPHTaitwaterGPU::init_style()
if (atom->molecular != Atom::ATOMIC) maxspecial = atom->maxspecial;
int mnf = 5e-2 * neighbor->oneatom;
int success =
sph_taitwater_gpu_init(atom->ntypes + 1, cutsq, cut, viscosity,
sph_taitwater_gpu_init(atom->ntypes + 1, cutsq, cut, viscosity, atom->mass,
rho0, soundspeed, B, domain->dimension, force->special_lj,
atom->nlocal, atom->nlocal + atom->nghost,
mnf, maxspecial, cell_size, gpu_mode, screen);