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

Fixed bugs with sph gpu pair styles

Browse Source
This commit is contained in:
Trung Nguyen
2024-03-18 13:27:21 -05:00
parent 207a14f351
commit 6a28e8d5f6
6 changed files with 112 additions and 99 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
lib/gpu/lal_sph_heatconduction.cu
View File

@ -29,23 +29,23 @@ _texture_2d( vel_tex,int4);
#if (SHUFFLE_AVAIL == 0)
#define store_dE(dEacc, ii, inum, tid, t_per_atom, offset, dE) \
#define store_dE(dEacc, ii, inum, tid, t_per_atom, offset, i, dE) \
if (t_per_atom>1) { \
simdsync(); \
simd_reduce_add1(t_per_atom, red_acc, offset, tid, dEacc); \
} \
if (offset==0 && ii<inum) { \
dE[ii]=dEacc; \
dE[i]=dEacc; \
}
#else
#define store_drhoE(dEacc, ii, inum, tid, t_per_atom, offset, dE) \
#define store_drhoE(dEacc, ii, inum, tid, t_per_atom, offset, i, dE) \
if (t_per_atom>1) { \
for (unsigned int s=t_per_atom/2; s>0; s>>=1) { \
dEacc += shfl_down(dEacc, s, t_per_atom); \
} \
} \
if (offset==0 && ii<inum) { \
dE[ii]=dEacc; \
dE[i]=dEacc; \
}
#endif
@ -66,7 +66,7 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int dimension, const int t_per_atom) {
int tid, ii, offset;
int tid, ii, offset, i;
atom_info(t_per_atom,ii,tid,offset);
int n_stride;
@ -77,7 +77,7 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
acctyp dEacc = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -140,7 +140,7 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
} // for nbor
} // if ii
store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,dE);
store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,i,dE);
}
__kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
@ -157,7 +157,7 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int dimension, const int t_per_atom) {
int tid, ii, offset;
int tid, ii, offset, i;
atom_info(t_per_atom,ii,tid,offset);
#ifndef ONETYPE
@ -180,7 +180,7 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
acctyp dEacc = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -204,7 +204,7 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
#endif
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
int jtype = jx.w;
int jtype=jx.w;
#ifndef ONETYPE
int mtype=itype+jx.w;
const numtyp cutsq_p=coeff[mtype].z;
@ -252,6 +252,6 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
} // for nbor
} // if ii
store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,dE);
store_drhoE(dEacc,ii,inum,tid,t_per_atom,offset,i,dE);
}

45
lib/gpu/lal_sph_lj.cu
View File

@ -29,17 +29,18 @@ _texture_2d( vel_tex,int4);
#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; \
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, 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[i]=drhoEacc.x; \
drhoE[i+inum]=drhoEacc.y; \
}
#else
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE) \
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, 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); \
@ -47,7 +48,8 @@ _texture_2d( vel_tex,int4);
} \
} \
if (offset==0 && ii<inum) { \
drhoE[ii]=drhoEacc; \
drhoE[i]=drhoEacc.x; \
drhoE[i+inum]=drhoEacc.y; \
}
#endif
@ -105,12 +107,12 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp2 *restrict drhoE,
__global acctyp *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int dimension, const int t_per_atom) {
int tid, ii, offset;
int tid, ii, offset, i;
atom_info(t_per_atom,ii,tid,offset);
int n_stride;
@ -124,10 +126,10 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
for (int i=0; i<6; i++) virial[i]=(acctyp)0;
}
acctyp2 drhoEacc;
drhoEacc.x = drhoEacc.x = (acctyp)0;
drhoEacc.x = drhoEacc.y = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -246,7 +248,7 @@ __kernel void k_sph_lj(const __global numtyp4 *restrict x_,
} // 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);
store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
}
__kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
@ -258,12 +260,12 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp2 *restrict drhoE,
__global acctyp *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int dimension, const int t_per_atom) {
int tid, ii, offset;
int tid, ii, offset, i;
atom_info(t_per_atom,ii,tid,offset);
#ifndef ONETYPE
@ -289,10 +291,10 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
for (int i=0; i<6; i++) virial[i]=(acctyp)0;
}
acctyp2 drhoEacc;
drhoEacc.x = drhoEacc.x = (acctyp)0;
drhoEacc.x = drhoEacc.y = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -325,7 +327,7 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
#endif
numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
int jtype = jx.w;
int jtype=jx.w;
#ifndef ONETYPE
int mtype=itype+jx.w;
const numtyp cutsq_p=coeff[mtype].z; // cutsq[itype][jtype];
@ -415,12 +417,11 @@ __kernel void k_sph_lj_fast(const __global numtyp4 *restrict x_,
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);
store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
}

64
lib/gpu/lal_sph_taitwater.cu
View File

@ -29,17 +29,18 @@ _texture_2d( vel_tex,int4);
#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; \
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, 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[i]=drhoEacc.x; \
drhoE[i+inum]=drhoEacc.y; \
}
#else
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, drhoE) \
#define store_drhoE(drhoEacc, ii, inum, tid, t_per_atom, offset, i, 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); \
@ -47,7 +48,8 @@ _texture_2d( vel_tex,int4);
} \
} \
if (offset==0 && ii<inum) { \
drhoE[ii]=drhoEacc; \
drhoE[i]=drhoEacc.x; \
drhoE[i+inum]=drhoEacc.y; \
}
#endif
@ -61,12 +63,12 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp2 *restrict drhoE,
__global acctyp *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int dimension, const int t_per_atom) {
int tid, ii, offset;
int tid, ii, offset, i;
atom_info(t_per_atom,ii,tid,offset);
int n_stride;
@ -80,10 +82,10 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
for (int i=0; i<6; i++) virial[i]=(acctyp)0;
}
acctyp2 drhoEacc;
drhoEacc.x = drhoEacc.x = (acctyp)0;
drhoEacc.x = drhoEacc.y = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -163,9 +165,8 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
numtyp fvisc = (numtyp)0;
if (delVdotDelR < (numtyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numtyp)0.01 * h * h);
fvisc = -coeffx * (soundspeed_itype
+ soundspeed_jtype) * mu / (rhoi + rhoj);
}
fvisc = -coeffx * (soundspeed_itype + soundspeed_jtype) * mu / (rhoi + rhoj);
}
// total pair force & thermal energy increment
numtyp force = -mass_itype * mass_jtype * (fi + fj + fvisc) * wfd;
@ -176,15 +177,11 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
f.z+=delz*force;
// and change in density, drho[i]
drhoEacc.x += mass_jtype* delVdotDelR * wfd;
drhoEacc.x += mass_jtype * 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;
@ -198,7 +195,7 @@ __kernel void k_sph_taitwater(const __global numtyp4 *restrict x_,
} // 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);
store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
}
__kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
@ -210,12 +207,12 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
const __global int * dev_packed,
__global acctyp3 *restrict ans,
__global acctyp *restrict engv,
__global acctyp2 *restrict drhoE,
__global acctyp *restrict drhoE,
const int eflag, const int vflag,
const int inum, const int nbor_pitch,
const __global numtyp4 *restrict v_,
const int dimension, const int t_per_atom) {
int tid, ii, offset;
int tid, ii, offset, i;
atom_info(t_per_atom,ii,tid,offset);
#ifndef ONETYPE
@ -245,10 +242,10 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
for (int i=0; i<6; i++) virial[i]=(acctyp)0;
}
acctyp2 drhoEacc;
drhoEacc.x = drhoEacc.x = (acctyp)0;
drhoEacc.x = drhoEacc.y = (acctyp)0;
if (ii<inum) {
int i, numj, nbor, nbor_end;
int numj, nbor, nbor_end;
nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
n_stride,nbor_end,nbor);
@ -337,8 +334,7 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
numtyp fvisc = (numtyp)0;
if (delVdotDelR < (numtyp)0) {
numtyp mu = h * delVdotDelR / (rsq + (numtyp)0.01 * h * h);
fvisc = -coeffx * (soundspeed_itype
+ soundspeed_jtype) * mu / (rhoi + rhoj);
fvisc = -coeffx * (soundspeed_itype + soundspeed_jtype) * mu / (rhoi + rhoj);
}
// total pair force & thermal energy increment
@ -349,16 +345,12 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
f.y+=dely*force;
f.z+=delz*force;
// and change in density
// and change in density, drho[i]
drhoEacc.x += mass_jtype * delVdotDelR * wfd;
// change in thermal energy
// 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;
@ -372,6 +364,6 @@ __kernel void k_sph_taitwater_fast(const __global numtyp4 *restrict x_,
} // 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);
store_drhoE(drhoEacc,ii,inum,tid,t_per_atom,offset,i,drhoE);
}

4
src/GPU/pair_sph_heatconduction_gpu.cpp
View File

@ -113,7 +113,7 @@ void PairSPHHeatConductionGPU::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);
cpu_time, success, atom->vest);
} else {
inum = list->inum;
ilist = list->ilist;
@ -122,7 +122,7 @@ void PairSPHHeatConductionGPU::compute(int eflag, int vflag)
sph_heatconduction_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->tag, atom->vest);
}
if (!success) error->one(FLERR, "Insufficient memory on accelerator");

28
src/GPU/pair_sph_lj_gpu.cpp
View File

@ -114,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);
cpu_time, success, atom->vest);
} else {
inum = list->inum;
ilist = list->ilist;
@ -123,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->tag, atom->vest);
}
if (!success) error->one(FLERR, "Insufficient memory on accelerator");
@ -136,21 +136,21 @@ void PairSPHLJGPU::compute(int eflag, int vflag)
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;
}
for (int i = 0; i < nlocal; i++)
drho[i] += drhoE_ptr[i];
drhoE_ptr += nlocal;
for (int i = 0; i < nlocal; i++)
desph[i] += drhoE_ptr[i];
} 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;
}
for (int i = 0; i < nlocal; i++)
drho[i] += drhoE_ptr[i];
drhoE_ptr += nlocal;
for (int i = 0; i < nlocal; i++)
desph[i] += drhoE_ptr[i];
}
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)

48
src/GPU/pair_sph_taitwater_gpu.cpp
View File

@ -18,6 +18,7 @@
#include "pair_sph_taitwater_gpu.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
@ -85,6 +86,25 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
{
ev_init(eflag, vflag);
// check consistency of pair coefficients
if (first) {
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = 1; i <= atom->ntypes; i++) {
if (cutsq[i][j] > 1.e-32) {
if (!setflag[i][i] || !setflag[j][j]) {
if (comm->me == 0) {
printf(
"SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n",
i, j, sqrt(cutsq[i][j]));
}
}
}
}
}
first = 0;
}
int nall = atom->nlocal + atom->nghost;
int inum, host_start;
@ -110,7 +130,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);
cpu_time, success, atom->vest);
} else {
inum = list->inum;
ilist = list->ilist;
@ -118,7 +138,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->tag, atom->vest);
}
if (!success) error->one(FLERR, "Insufficient memory on accelerator");
@ -131,21 +151,21 @@ void PairSPHTaitwaterGPU::compute(int eflag, int vflag)
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;
}
for (int i = 0; i < nlocal; i++)
drho[i] += drhoE_ptr[i];
drhoE_ptr += nlocal;
for (int i = 0; i < nlocal; i++)
desph[i] += drhoE_ptr[i];
} 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;
}
for (int i = 0; i < nlocal; i++)
drho[i] += drhoE_ptr[i];
drhoE_ptr += nlocal;
for (int i = 0; i < nlocal; i++)
desph[i] += drhoE_ptr[i];
}
if (atom->molecular != Atom::ATOMIC && neighbor->ago == 0)