whitespace

This commit is contained in:
Stan Moore
2024-10-24 09:43:36 -06:00
parent a844959a07
commit f18850397d
13 changed files with 335 additions and 341 deletions

View File

@ -56,6 +56,7 @@ FixCMAPKokkos<DeviceType>::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) :
datamask_modify = F_MASK; datamask_modify = F_MASK;
// allocate memory for CMAP data // allocate memory for CMAP data
memoryKK->create_kokkos(k_g_axis,g_axis,CMAPDIM,"cmap:g_axis"); memoryKK->create_kokkos(k_g_axis,g_axis,CMAPDIM,"cmap:g_axis");
memoryKK->create_kokkos(k_cmapgrid,cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:grid"); memoryKK->create_kokkos(k_cmapgrid,cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:grid");
memoryKK->create_kokkos(k_d1cmapgrid,d1cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:d1grid"); memoryKK->create_kokkos(k_d1cmapgrid,d1cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:d1grid");
@ -69,6 +70,7 @@ FixCMAPKokkos<DeviceType>::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) :
d_d12cmapgrid = k_d12cmapgrid.template view<DeviceType>(); d_d12cmapgrid = k_d12cmapgrid.template view<DeviceType>();
// read and setup CMAP data // read and setup CMAP data
read_grid_map(arg[3]); read_grid_map(arg[3]);
int i = 0; int i = 0;
@ -88,6 +90,7 @@ FixCMAPKokkos<DeviceType>::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) :
for( int i=0 ; i<CMAPMAX ; i++ ) { for( int i=0 ; i<CMAPMAX ; i++ ) {
// pre-compute the derivatives of the maps // pre-compute the derivatives of the maps
set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]); set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]);
for( int j=0 ; j<CMAPDIM ; j++ ) { for( int j=0 ; j<CMAPDIM ; j++ ) {
@ -111,7 +114,6 @@ FixCMAPKokkos<DeviceType>::FixCMAPKokkos(LAMMPS *lmp, int narg, char **arg) :
k_d1cmapgrid.template sync<DeviceType>(); k_d1cmapgrid.template sync<DeviceType>();
k_d2cmapgrid.template sync<DeviceType>(); k_d2cmapgrid.template sync<DeviceType>();
k_d12cmapgrid.template sync<DeviceType>(); k_d12cmapgrid.template sync<DeviceType>();
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -136,7 +138,6 @@ FixCMAPKokkos<DeviceType>::~FixCMAPKokkos()
memoryKK->destroy_kokkos(k_crossterm_atom5,crossterm_atom5); memoryKK->destroy_kokkos(k_crossterm_atom5,crossterm_atom5);
memoryKK->destroy_kokkos(d_crosstermlist); memoryKK->destroy_kokkos(d_crosstermlist);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -148,6 +149,7 @@ void FixCMAPKokkos<DeviceType>::init()
error->all(FLERR,"Cannot yet use respa with Kokkos"); error->all(FLERR,"Cannot yet use respa with Kokkos");
// on KOKKOS, allocate enough for all crossterms on each GPU to avoid grow operation in device code // on KOKKOS, allocate enough for all crossterms on each GPU to avoid grow operation in device code
maxcrossterm = ncmap; maxcrossterm = ncmap;
memoryKK->create_kokkos(d_crosstermlist,maxcrossterm,CMAPMAX,"cmap:crosstermlist"); memoryKK->create_kokkos(d_crosstermlist,maxcrossterm,CMAPMAX,"cmap:crosstermlist");
} }
@ -159,7 +161,6 @@ void FixCMAPKokkos<DeviceType>::init()
template<class DeviceType> template<class DeviceType>
void FixCMAPKokkos<DeviceType>::pre_neighbor() void FixCMAPKokkos<DeviceType>::pre_neighbor()
{ {
atomKK->sync(execution_space,X_MASK); atomKK->sync(execution_space,X_MASK);
d_x = atomKK->k_x.view<DeviceType>(); d_x = atomKK->k_x.view<DeviceType>();
int nlocal = atomKK->nlocal; int nlocal = atomKK->nlocal;
@ -179,14 +180,12 @@ void FixCMAPKokkos<DeviceType>::pre_neighbor()
copymode = 1; copymode = 1;
Kokkos::parallel_scan(Kokkos::RangePolicy<DeviceType,TagFixCmapPreNeighbor>(0,nlocal),*this,ncrosstermlist); Kokkos::parallel_scan(Kokkos::RangePolicy<DeviceType,TagFixCmapPreNeighbor>(0,nlocal),*this,ncrosstermlist);
copymode = 0; copymode = 0;
} }
template<class DeviceType> template<class DeviceType>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
void FixCMAPKokkos<DeviceType>::operator()(TagFixCmapPreNeighbor, const int i, int &l_ncrosstermlist, const bool is_final ) const void FixCMAPKokkos<DeviceType>::operator()(TagFixCmapPreNeighbor, const int i, int &l_ncrosstermlist, const bool is_final ) const
{ {
for( int m = 0; m < d_num_crossterm(i); m++) { for( int m = 0; m < d_num_crossterm(i); m++) {
int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_crossterm_atom1(i,m),map_style,k_map_array,k_map_hash); int atom1 = AtomKokkos::map_kokkos<DeviceType>(d_crossterm_atom1(i,m),map_style,k_map_array,k_map_hash);
@ -297,6 +296,7 @@ void FixCMAPKokkos<DeviceType>::operator()(TagFixCmapPostForce, const int n, dou
double vb45z = d_x(i4,2) - d_x(i5,2); double vb45z = d_x(i4,2) - d_x(i5,2);
// calculate normal vectors for planes that define the dihedral angles // calculate normal vectors for planes that define the dihedral angles
double a1x = vb12y*vb23z - vb12z*vb23y; double a1x = vb12y*vb23z - vb12z*vb23y;
double a1y = vb12z*vb23x - vb12x*vb23z; double a1y = vb12z*vb23x - vb12x*vb23z;
double a1z = vb12x*vb23y - vb12y*vb23x; double a1z = vb12x*vb23y - vb12y*vb23x;
@ -325,6 +325,7 @@ void FixCMAPKokkos<DeviceType>::operator()(TagFixCmapPostForce, const int n, dou
if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) return; if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) return;
// vectors needed to calculate the cross-term dihedral angles // vectors needed to calculate the cross-term dihedral angles
double dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z; double dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z;
double dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z; double dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z;
double dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z; double dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z;
@ -388,8 +389,8 @@ void FixCMAPKokkos<DeviceType>::operator()(TagFixCmapPostForce, const int n, dou
bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs,E,dEdPhi,dEdPsi); bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs,E,dEdPhi,dEdPsi);
// sum up cmap energy contributions // sum up cmap energy contributions
// needed for compute_scalar() // needed for compute_scalar()
double engfraction = 0.2 * E; double engfraction = 0.2 * E;
if (i1 < nlocal) ecmapKK += engfraction; if (i1 < nlocal) ecmapKK += engfraction;
if (i2 < nlocal) ecmapKK += engfraction; if (i2 < nlocal) ecmapKK += engfraction;
@ -479,6 +480,7 @@ void FixCMAPKokkos<DeviceType>::grow_arrays(int nmax)
k_crossterm_atom5.template sync<LMPHostType>(); k_crossterm_atom5.template sync<LMPHostType>();
// force reallocation on host // force reallocation on host
k_num_crossterm.template modify<LMPHostType>(); k_num_crossterm.template modify<LMPHostType>();
k_crossterm_type.template modify<LMPHostType>(); k_crossterm_type.template modify<LMPHostType>();
k_crossterm_atom1.template modify<LMPHostType>(); k_crossterm_atom1.template modify<LMPHostType>();
@ -877,8 +879,6 @@ void FixCMAPKokkos<DeviceType>::bc_interpol(double x1, double x2, int low1, int
dEdPsi *= (180.0/MY_PI/CMAPDX); dEdPsi *= (180.0/MY_PI/CMAPDX);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
return local index of atom J or any of its images that is closest to atom I return local index of atom J or any of its images that is closest to atom I
if J is not a valid index like -1, just return it if J is not a valid index like -1, just return it
@ -917,7 +917,6 @@ int FixCMAPKokkos<DeviceType>::closest_image(const int i, int j) const
return closest; return closest;
} }
namespace LAMMPS_NS { namespace LAMMPS_NS {
template class FixCMAPKokkos<LMPDeviceType>; template class FixCMAPKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU #ifdef LMP_KOKKOS_GPU

View File

@ -120,7 +120,6 @@ void FixNVELimitKokkos<DeviceType>::initial_integrate(int /*vflag*/)
ncount += d_ncount; ncount += d_ncount;
atomKK->modified(execution_space, X_MASK | V_MASK ); atomKK->modified(execution_space, X_MASK | V_MASK );
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -190,7 +189,6 @@ void FixNVELimitKokkos<DeviceType>::final_integrate()
ncount += d_ncount; ncount += d_ncount;
atomKK->modified(execution_space, V_MASK ); atomKK->modified(execution_space, V_MASK );
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

View File

@ -50,7 +50,6 @@ FixRecenterKokkos<DeviceType>::FixRecenterKokkos(LAMMPS *lmp, int narg, char **a
template<class DeviceType> template<class DeviceType>
void FixRecenterKokkos<DeviceType>::initial_integrate(int /*vflag*/) void FixRecenterKokkos<DeviceType>::initial_integrate(int /*vflag*/)
{ {
atomKK->sync(execution_space,datamask_read); atomKK->sync(execution_space,datamask_read);
int nlocal = atomKK->nlocal; int nlocal = atomKK->nlocal;
if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
@ -121,7 +120,6 @@ void FixRecenterKokkos<DeviceType>::initial_integrate(int /*vflag*/)
atomKK->modified(execution_space,datamask_modify); atomKK->modified(execution_space,datamask_modify);
} }
namespace LAMMPS_NS { namespace LAMMPS_NS {
template class FixRecenterKokkos<LMPDeviceType>; template class FixRecenterKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU #ifdef LMP_KOKKOS_GPU

View File

@ -91,8 +91,8 @@ void FixWallRegionKokkos<DeviceType>::post_force(int vflag)
// initilize ewall after region->prematch(), // initilize ewall after region->prematch(),
// so a dynamic region can access last timestep values // so a dynamic region can access last timestep values
// energy intialize. // energy intialize
// eflag is used to track whether wall energies have been communicated. // eflag is used to track whether wall energies have been communicated
eflag = 0; eflag = 0;
double result[10]; double result[10];
@ -330,7 +330,6 @@ template <class DeviceType>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
void FixWallRegionKokkos<DeviceType>::v_tally(value_type result, int i, double *v) const void FixWallRegionKokkos<DeviceType>::v_tally(value_type result, int i, double *v) const
{ {
if (vflag_global) { if (vflag_global) {
result[4] += v[0]; result[4] += v[0];
result[5] += v[1]; result[5] += v[1];
@ -348,7 +347,6 @@ void FixWallRegionKokkos<DeviceType>::v_tally(value_type result, int i, double *
Kokkos::atomic_add(&(d_vatom(i,4)),v[4]); Kokkos::atomic_add(&(d_vatom(i,4)),v[4]);
Kokkos::atomic_add(&(d_vatom(i,5)),v[5]); Kokkos::atomic_add(&(d_vatom(i,5)),v[5]);
} }
} }
namespace LAMMPS_NS { namespace LAMMPS_NS {

View File

@ -71,7 +71,6 @@ double GroupKokkos<DeviceType>::mass(int igroup)
return all; return all;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
compute the center-of-mass coords of group of atoms compute the center-of-mass coords of group of atoms
masstotal = total mass masstotal = total mass

View File

@ -23,7 +23,7 @@ template<class DeviceType>
class GroupKokkos : public Group { class GroupKokkos : public Group {
public: public:
GroupKokkos(class LAMMPS *); GroupKokkos(class LAMMPS *);
double mass(int); // total mass of atoms in group double mass(int); // total mass of atoms in group
void xcm(int, double, double *); // center-of-mass coords of group void xcm(int, double, double *); // center-of-mass coords of group
}; };

View File

@ -271,7 +271,6 @@ namespace MathSpecialKokkos {
return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
} }
} // namespace MathSpecialKokkos } // namespace MathSpecialKokkos
} // namespace LAMMPS_NS } // namespace LAMMPS_NS

View File

@ -61,43 +61,44 @@ class RegBlockKokkos : public RegBlock, public KokkosBase {
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
int surface_kokkos(double x, double y, double z, double cutoff) int surface_kokkos(double x, double y, double z, double cutoff)
{ {
int ncontact; int ncontact;
double xs, ys, zs; double xs, ys, zs;
double xnear[3], xorig[3]; double xnear[3], xorig[3];
if (dynamic) { if (dynamic) {
xorig[0] = x; xorig[1] = y; xorig[2] = z; xorig[0] = x; xorig[1] = y; xorig[2] = z;
inverse_transform(x, y, z); inverse_transform(x, y, z);
}
xnear[0] = x; xnear[1] = y; xnear[2] = z;
if (!openflag) {
if (interior)
ncontact = surface_interior_kokkos(xnear, cutoff);
else
ncontact = surface_exterior_kokkos(xnear, cutoff);
} else {
// one of surface_int/ext() will return 0
// so no need to worry about offset of contact indices
ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff);
}
if (rotateflag && ncontact) {
for (int i = 0; i < ncontact; i++) {
xs = xnear[0] - d_contact[i].delx;
ys = xnear[1] - d_contact[i].dely;
zs = xnear[2] - d_contact[i].delz;
forward_transform(xs, ys, zs);
d_contact[i].delx = xorig[0] - xs;
d_contact[i].dely = xorig[1] - ys;
d_contact[i].delz = xorig[2] - zs;
} }
xnear[0] = x; xnear[1] = y; xnear[2] = z;
if (!openflag) {
if (interior)
ncontact = surface_interior_kokkos(xnear, cutoff);
else
ncontact = surface_exterior_kokkos(xnear, cutoff);
} else {
// one of surface_int/ext() will return 0
// so no need to worry about offset of contact indices
ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff);
}
if (rotateflag && ncontact) {
for (int i = 0; i < ncontact; i++) {
xs = xnear[0] - d_contact[i].delx;
ys = xnear[1] - d_contact[i].dely;
zs = xnear[2] - d_contact[i].delz;
forward_transform(xs, ys, zs);
d_contact[i].delx = xorig[0] - xs;
d_contact[i].dely = xorig[1] - ys;
d_contact[i].delz = xorig[2] - zs;
}
}
return ncontact;
} }
return ncontact;
}
Kokkos::View<Contact*, DeviceType> d_contact; Kokkos::View<Contact*, DeviceType> d_contact;
private: private:
@ -106,317 +107,316 @@ class RegBlockKokkos : public RegBlock, public KokkosBase {
typename AT::t_x_array_randomread d_x; typename AT::t_x_array_randomread d_x;
typename AT::t_int_1d_randomread d_mask; typename AT::t_int_1d_randomread d_mask;
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
int surface_interior_kokkos(double *x, double cutoff) int surface_interior_kokkos(double *x, double cutoff)
{ {
double delta; double delta;
// x is exterior to block // x is exterior to block
if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0; if (x[0] < xlo || x[0] > xhi || x[1] < ylo || x[1] > yhi || x[2] < zlo || x[2] > zhi) return 0;
// x is interior to block or on its surface // x is interior to block or on its surface
int n = 0; int n = 0;
delta = x[0] - xlo; delta = x[0] - xlo;
if (delta < cutoff && !open_faces[0]) { if (delta < cutoff && !open_faces[0]) {
d_contact[n].r = delta; d_contact[n].r = delta;
d_contact[n].delx = delta; d_contact[n].delx = delta;
d_contact[n].dely = d_contact[n].delz = 0.0; d_contact[n].dely = d_contact[n].delz = 0.0;
d_contact[n].radius = 0; d_contact[n].radius = 0;
d_contact[n].iwall = 0; d_contact[n].iwall = 0;
n++; n++;
} }
delta = xhi - x[0]; delta = xhi - x[0];
if (delta < cutoff && !open_faces[1]) { if (delta < cutoff && !open_faces[1]) {
d_contact[n].r = delta; d_contact[n].r = delta;
d_contact[n].delx = -delta; d_contact[n].delx = -delta;
d_contact[n].dely = d_contact[n].delz = 0.0; d_contact[n].dely = d_contact[n].delz = 0.0;
d_contact[n].radius = 0; d_contact[n].radius = 0;
d_contact[n].iwall = 1; d_contact[n].iwall = 1;
n++; n++;
}
delta = x[1] - ylo;
if (delta < cutoff && !open_faces[2]) {
d_contact[n].r = delta;
d_contact[n].dely = delta;
d_contact[n].delx = d_contact[n].delz = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 2;
n++;
}
delta = yhi - x[1];
if (delta < cutoff && !open_faces[3]) {
d_contact[n].r = delta;
d_contact[n].dely = -delta;
d_contact[n].delx = d_contact[n].delz = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 3;
n++;
}
delta = x[2] - zlo;
if (delta < cutoff && !open_faces[4]) {
d_contact[n].r = delta;
d_contact[n].delz = delta;
d_contact[n].delx = d_contact[n].dely = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 4;
n++;
}
delta = zhi - x[2];
if (delta < cutoff && !open_faces[5]) {
d_contact[n].r = delta;
d_contact[n].delz = -delta;
d_contact[n].delx = d_contact[n].dely = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 5;
n++;
}
return n;
} }
delta = x[1] - ylo; KOKKOS_INLINE_FUNCTION
if (delta < cutoff && !open_faces[2]) { int surface_exterior_kokkos(double *x, double cutoff)
d_contact[n].r = delta; {
d_contact[n].dely = delta; double xp, yp, zp;
d_contact[n].delx = d_contact[n].delz = 0.0; double xc, yc, zc, dist, mindist;
d_contact[n].radius = 0;
d_contact[n].iwall = 2;
n++;
}
delta = yhi - x[1];
if (delta < cutoff && !open_faces[3]) {
d_contact[n].r = delta;
d_contact[n].dely = -delta;
d_contact[n].delx = d_contact[n].delz = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 3;
n++;
}
delta = x[2] - zlo; // x is far enough from block that there is no contact
if (delta < cutoff && !open_faces[4]) { // x is interior to block
d_contact[n].r = delta;
d_contact[n].delz = delta;
d_contact[n].delx = d_contact[n].dely = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 4;
n++;
}
delta = zhi - x[2];
if (delta < cutoff && !open_faces[5]) {
d_contact[n].r = delta;
d_contact[n].delz = -delta;
d_contact[n].delx = d_contact[n].dely = 0.0;
d_contact[n].radius = 0;
d_contact[n].iwall = 5;
n++;
}
return n; if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff ||
} x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff)
return 0;
if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0;
KOKKOS_INLINE_FUNCTION // x is exterior to block or on its surface
int surface_exterior_kokkos(double *x, double cutoff) // xp,yp,zp = point on surface of block that x is closest to
{ // could be edge or corner pt of block
double xp, yp, zp; // do not add contact point if r >= cutoff
double xc, yc, zc, dist, mindist;
// x is far enough from block that there is no contact if (!openflag) {
// x is interior to block if (x[0] < xlo)
xp = xlo;
if (x[0] <= xlo - cutoff || x[0] >= xhi + cutoff || x[1] <= ylo - cutoff || else if (x[0] > xhi)
x[1] >= yhi + cutoff || x[2] <= zlo - cutoff || x[2] >= zhi + cutoff) xp = xhi;
return 0; else
if (x[0] > xlo && x[0] < xhi && x[1] > ylo && x[1] < yhi && x[2] > zlo && x[2] < zhi) return 0; xp = x[0];
if (x[1] < ylo)
// x is exterior to block or on its surface yp = ylo;
// xp,yp,zp = point on surface of block that x is closest to else if (x[1] > yhi)
// could be edge or corner pt of block yp = yhi;
// do not add contact point if r >= cutoff else
yp = x[1];
if (!openflag) { if (x[2] < zlo)
if (x[0] < xlo) zp = zlo;
xp = xlo; else if (x[2] > zhi)
else if (x[0] > xhi) zp = zhi;
xp = xhi; else
else zp = x[2];
xp = x[0]; } else {
if (x[1] < ylo) mindist = MAXDOUBLEINT;
yp = ylo; for (int i = 0; i < 6; i++) {
else if (x[1] > yhi) if (open_faces[i]) continue;
yp = yhi; dist = find_closest_point(i, x, xc, yc, zc);
else if (dist < mindist) {
yp = x[1]; xp = xc;
if (x[2] < zlo) yp = yc;
zp = zlo; zp = zc;
else if (x[2] > zhi) mindist = dist;
zp = zhi; }
else
zp = x[2];
} else {
mindist = MAXDOUBLEINT;
for (int i = 0; i < 6; i++) {
if (open_faces[i]) continue;
dist = find_closest_point(i, x, xc, yc, zc);
if (dist < mindist) {
xp = xc;
yp = yc;
zp = zc;
mindist = dist;
} }
} }
add_contact(0, x, xp, yp, zp);
d_contact[0].iwall = 0;
if (d_contact[0].r < cutoff) return 1;
return 0;
} }
add_contact(0, x, xp, yp, zp); KOKKOS_INLINE_FUNCTION
d_contact[0].iwall = 0; void add_contact(int n, double *x, double xp, double yp, double zp)
if (d_contact[0].r < cutoff) return 1; {
return 0; double delx = x[0] - xp;
} double dely = x[1] - yp;
double delz = x[2] - zp;
KOKKOS_INLINE_FUNCTION d_contact[n].r = sqrt(delx * delx + dely * dely + delz * delz);
void add_contact(int n, double *x, double xp, double yp, double zp) d_contact[n].radius = 0;
{ d_contact[n].delx = delx;
double delx = x[0] - xp; d_contact[n].dely = dely;
double dely = x[1] - yp; d_contact[n].delz = delz;
double delz = x[2] - zp;
d_contact[n].r = sqrt(delx * delx + dely * dely + delz * delz);
d_contact[n].radius = 0;
d_contact[n].delx = delx;
d_contact[n].dely = dely;
d_contact[n].delz = delz;
}
KOKKOS_INLINE_FUNCTION
int k_inside(double x, double y, double z) const
{
if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi)
return 1;
return 0;
}
KOKKOS_INLINE_FUNCTION
void forward_transform(double &x, double &y, double &z) const
{
if (rotateflag) rotate(x, y, z, theta);
if (moveflag) {
x += dx;
y += dy;
z += dz;
} }
}
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
void inverse_transform(double &x, double &y, double &z) const int k_inside(double x, double y, double z) const
{ {
if (moveflag) { if (x >= xlo && x <= xhi && y >= ylo && y <= yhi && z >= zlo && z <= zhi)
x -= dx;
y -= dy;
z -= dz;
}
if (rotateflag) rotate(x,y,z,-theta);
}
KOKKOS_INLINE_FUNCTION
void rotate(double &x, double &y, double &z, double angle) const
{
double a[3],b[3],c[3],d[3],disp[3];
double sine = sin(angle);
double cosine = cos(angle);
d[0] = x - point[0];
d[1] = y - point[1];
d[2] = z - point[2];
double x0dotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
c[0] = x0dotr * runit[0];
c[1] = x0dotr * runit[1];
c[2] = x0dotr * runit[2];
a[0] = d[0] - c[0];
a[1] = d[1] - c[1];
a[2] = d[2] - c[2];
b[0] = runit[1]*a[2] - runit[2]*a[1];
b[1] = runit[2]*a[0] - runit[0]*a[2];
b[2] = runit[0]*a[1] - runit[1]*a[0];
disp[0] = a[0]*cosine + b[0]*sine;
disp[1] = a[1]*cosine + b[1]*sine;
disp[2] = a[2]*cosine + b[2]*sine;
x = point[0] + c[0] + disp[0];
y = point[1] + c[1] + disp[1];
z = point[2] + c[2] + disp[2];
}
KOKKOS_INLINE_FUNCTION
void point_on_line_segment(double *a, double *b, double *c, double *d)
{
double ba[3], ca[3];
sub3(b, a, ba);
sub3(c, a, ca);
double t = dot3(ca, ba) / dot3(ba, ba);
if (t <= 0.0) {
d[0] = a[0];
d[1] = a[1];
d[2] = a[2];
} else if (t >= 1.0) {
d[0] = b[0];
d[1] = b[1];
d[2] = b[2];
} else {
d[0] = a[0] + t * ba[0];
d[1] = a[1] + t * ba[1];
d[2] = a[2] + t * ba[2];
}
}
KOKKOS_INLINE_FUNCTION
double inside_face(double *xproj, int iface)
{
if (iface < 2) {
if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1;
} else if (iface < 4) {
if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo)))
return 1; return 1;
} else { return 0;
if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1;
} }
return 0; KOKKOS_INLINE_FUNCTION
} void forward_transform(double &x, double &y, double &z) const
{
KOKKOS_INLINE_FUNCTION if (rotateflag) rotate(x, y, z, theta);
double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) if (moveflag) {
{ x += dx;
double dot, d2, d2min; y += dy;
double xr[3], xproj[3], p[3]; z += dz;
xr[0] = x[0] - corners[i][0][0];
xr[1] = x[1] - corners[i][0][1];
xr[2] = x[2] - corners[i][0][2];
dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2];
xproj[0] = xr[0] - dot * face[i][0];
xproj[1] = xr[1] - dot * face[i][1];
xproj[2] = xr[2] - dot * face[i][2];
d2min = MAXDOUBLEINT;
// check if point projects inside of face
if (inside_face(xproj, i)) {
d2 = d2min = dot * dot;
xc = xproj[0] + corners[i][0][0];
yc = xproj[1] + corners[i][0][1];
zc = xproj[2] + corners[i][0][2];
// check each edge
} else {
point_on_line_segment(corners[i][0], corners[i][1], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
point_on_line_segment(corners[i][1], corners[i][2], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
point_on_line_segment(corners[i][2], corners[i][3], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
point_on_line_segment(corners[i][3], corners[i][0], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
} }
} }
return d2min; KOKKOS_INLINE_FUNCTION
} void inverse_transform(double &x, double &y, double &z) const
{
if (moveflag) {
x -= dx;
y -= dy;
z -= dz;
}
if (rotateflag) rotate(x,y,z,-theta);
}
KOKKOS_INLINE_FUNCTION
void rotate(double &x, double &y, double &z, double angle) const
{
double a[3],b[3],c[3],d[3],disp[3];
double sine = sin(angle);
double cosine = cos(angle);
d[0] = x - point[0];
d[1] = y - point[1];
d[2] = z - point[2];
double x0dotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
c[0] = x0dotr * runit[0];
c[1] = x0dotr * runit[1];
c[2] = x0dotr * runit[2];
a[0] = d[0] - c[0];
a[1] = d[1] - c[1];
a[2] = d[2] - c[2];
b[0] = runit[1]*a[2] - runit[2]*a[1];
b[1] = runit[2]*a[0] - runit[0]*a[2];
b[2] = runit[0]*a[1] - runit[1]*a[0];
disp[0] = a[0]*cosine + b[0]*sine;
disp[1] = a[1]*cosine + b[1]*sine;
disp[2] = a[2]*cosine + b[2]*sine;
x = point[0] + c[0] + disp[0];
y = point[1] + c[1] + disp[1];
z = point[2] + c[2] + disp[2];
}
KOKKOS_INLINE_FUNCTION
void point_on_line_segment(double *a, double *b, double *c, double *d)
{
double ba[3], ca[3];
sub3(b, a, ba);
sub3(c, a, ca);
double t = dot3(ca, ba) / dot3(ba, ba);
if (t <= 0.0) {
d[0] = a[0];
d[1] = a[1];
d[2] = a[2];
} else if (t >= 1.0) {
d[0] = b[0];
d[1] = b[1];
d[2] = b[2];
} else {
d[0] = a[0] + t * ba[0];
d[1] = a[1] + t * ba[1];
d[2] = a[2] + t * ba[2];
}
}
KOKKOS_INLINE_FUNCTION
double inside_face(double *xproj, int iface)
{
if (iface < 2) {
if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1;
} else if (iface < 4) {
if (xproj[0] > 0 && (xproj[0] < (xhi - xlo)) && xproj[2] > 0 && (xproj[2] < (zhi - zlo)))
return 1;
} else {
if (xproj[0] > 0 && xproj[0] < (xhi - xlo) && xproj[1] > 0 && xproj[1] < (yhi - ylo)) return 1;
}
return 0;
}
KOKKOS_INLINE_FUNCTION
double find_closest_point(int i, double *x, double &xc, double &yc, double &zc)
{
double dot, d2, d2min;
double xr[3], xproj[3], p[3];
xr[0] = x[0] - corners[i][0][0];
xr[1] = x[1] - corners[i][0][1];
xr[2] = x[2] - corners[i][0][2];
dot = face[i][0] * xr[0] + face[i][1] * xr[1] + face[i][2] * xr[2];
xproj[0] = xr[0] - dot * face[i][0];
xproj[1] = xr[1] - dot * face[i][1];
xproj[2] = xr[2] - dot * face[i][2];
d2min = MAXDOUBLEINT;
// check if point projects inside of face
if (inside_face(xproj, i)) {
d2 = d2min = dot * dot;
xc = xproj[0] + corners[i][0][0];
yc = xproj[1] + corners[i][0][1];
zc = xproj[2] + corners[i][0][2];
// check each edge
} else {
point_on_line_segment(corners[i][0], corners[i][1], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
point_on_line_segment(corners[i][1], corners[i][2], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
point_on_line_segment(corners[i][2], corners[i][3], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
point_on_line_segment(corners[i][3], corners[i][0], x, p);
d2 = (p[0] - x[0]) * (p[0] - x[0]) + (p[1] - x[1]) * (p[1] - x[1]) +
(p[2] - x[2]) * (p[2] - x[2]);
if (d2 < d2min) {
d2min = d2;
xc = p[0];
yc = p[1];
zc = p[2];
}
}
return d2min;
}
}; };

View File

@ -75,7 +75,6 @@ void RegSphereKokkos<DeviceType>::operator()(TagRegSphereMatchAll, const int &i)
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
namespace LAMMPS_NS { namespace LAMMPS_NS {

View File

@ -131,6 +131,7 @@ FixCMAP::~FixCMAP()
if (copymode) return; if (copymode) return;
// unregister callbacks to this fix from Atom class // unregister callbacks to this fix from Atom class
atom->delete_callback(id,Atom::GROW); atom->delete_callback(id,Atom::GROW);
atom->delete_callback(id,Atom::RESTART); atom->delete_callback(id,Atom::RESTART);

View File

@ -50,6 +50,7 @@ Region::Region(LAMMPS *lmp, int /*narg*/, char **arg) :
Region::~Region() Region::~Region()
{ {
if (copymode) return; if (copymode) return;
delete[] id; delete[] id;
delete[] style; delete[] style;
delete[] xstr; delete[] xstr;

View File

@ -262,6 +262,7 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
RegBlock::~RegBlock() RegBlock::~RegBlock()
{ {
if (copymode) return; if (copymode) return;
delete[] xlostr; delete[] xlostr;
delete[] xhistr; delete[] xhistr;
delete[] ylostr; delete[] ylostr;

View File

@ -102,6 +102,7 @@ RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
RegSphere::~RegSphere() RegSphere::~RegSphere()
{ {
if (copymode) return; if (copymode) return;
delete[] xstr; delete[] xstr;
delete[] ystr; delete[] ystr;
delete[] zstr; delete[] zstr;