diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 102e16a711b..06bb7a9ad1f 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -117,7 +117,7 @@ parenthesis: g = GPU, i = INTEL, k = KOKKOS, o = OPENMP, t = OPT. * :doc:`gran/hertz/history (o) ` * :doc:`gran/hooke (o) ` * :doc:`gran/hooke/history (ko) ` - * :doc:`granular ` + * :doc:`granular (k) ` * :doc:`granular/superellipsoid ` * :doc:`gw (o) ` * :doc:`gw/zbl (o) ` diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index e2d20354108..d5cc90c956f 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -1,8 +1,11 @@ .. index:: pair_style granular +.. index:: pair_style granular/kk pair_style granular command =========================== +Accelerator Variants: *granular/kk* + Syntax """""" @@ -1054,6 +1057,10 @@ These extra quantities can be accessed by the :doc:`compute pair/local ---------- +.. include:: accel_styles.rst + +---------- + Restrictions """""""""""" @@ -1061,6 +1068,38 @@ This pair style is part of the GRANULAR package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. +The *granular/kk* style also requires the KOKKOS package. Its model +support matrix is: + +.. list-table:: + :header-rows: 1 + + * - Category + - Supported models + - Unsupported models + * - Normal + - *hooke*, *hertz*, *hertz/material* + - *mdr* and other normal models + * - Damping + - *none*, *velocity*, *mass_velocity*, *viscoelastic*, *tsuji*, + *coeff_restitution* + - *mdr* and other damping models + * - Tangential + - *none*, *linear_nohistory*, *linear_history*, *mindlin*, + *mindlin/force* + - *mindlin_rescale*, *mindlin_rescale/force*, and other tangential models + * - Rolling + - *none*, *sds* + - other rolling models + * - Twisting + - *none*, *sds*, *marshall* + - other twisting models + * - Heat + - *none*, *radius*, *area* + - other heat models + +Models requiring nondefault contact-history transfer are unsupported. + This pair style requires that atoms store per-particle radius, torque, and angular velocity (omega) as defined by the :doc:`atom_style sphere `. diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index bab3889de43..d330c26a7e3 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -249,10 +249,8 @@ void GranularModel::init() for (int i = 0; i < NSUBMODELS; i++) if (!sub_models[i]) construct_sub_model("none", (SubModelType) i); - // Must have valid normal, damping, and tangential models + // Must have a valid normal model if (normal_model->name == "none") error->all(FLERR, "Must specify normal granular model"); - if (damping_model->name == "none") error->all(FLERR, "Must specify damping granular model"); - if (tangential_model->name == "none") error->all(FLERR, "Must specify tangential granular model"); // Twisting, rolling, and heat are optional twisting_defined = rolling_defined = heat_defined = 1; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 32e13cad3f6..bcfc574dcfb 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -95,6 +95,8 @@ id_history = utils::strdup(std::string("NEIGH_HISTORY_GRANULAR") + std::to_strin PairGranular::~PairGranular() { + if (copymode) return; + delete[] svector; if (!fix_history) modify->delete_fix(id_dummy); @@ -492,7 +494,8 @@ void PairGranular::init_style() // this is so its order in the fix list is preserved if (use_history && fix_history == nullptr) { - fix_history = dynamic_cast(modify->replace_fix(id_dummy, fmt::format("{} all NEIGH_HISTORY {}", id_history, size_history),1)); + fix_history = dynamic_cast(modify->replace_fix( + id_dummy, fmt::format("{} all {} {}", id_history, history_fix_style(), size_history), 1)); fix_history->pair = this; } else if (use_history) { fix_history = dynamic_cast(modify->get_fix_by_id(id_history)); @@ -559,6 +562,13 @@ void PairGranular::init_style() MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world); } +/* ---------------------------------------------------------------------- */ + +const char *PairGranular::history_fix_style() const +{ + return "NEIGH_HISTORY"; +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 84a231a8377..6317eb88f10 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -72,8 +72,9 @@ class PairGranular : public Pair { double *mass_rigid; // rigid mass for owned+ghost atoms int nmax; // allocated size of mass_rigid - void allocate(); + virtual void allocate(); void prune_models(); + virtual const char *history_fix_style() const; int size_history; int heat_flag; diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 4e1429da7d7..e9fd963b1f7 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -416,6 +416,8 @@ action pair_gauss_kokkos.cpp pair_gauss.cpp action pair_gauss_kokkos.h pair_gauss.h action pair_gran_hooke_history_kokkos.cpp pair_gran_hooke_history.cpp action pair_gran_hooke_history_kokkos.h pair_gran_hooke_history.h +action pair_granular_kokkos.cpp pair_granular.cpp +action pair_granular_kokkos.h pair_granular.h action pair_gauss_cut_kokkos.cpp pair_gauss_cut.cpp action pair_gauss_cut_kokkos.h pair_gauss_cut.h action pair_hybrid_kokkos.cpp diff --git a/src/KOKKOS/pair_granular_kokkos.cpp b/src/KOKKOS/pair_granular_kokkos.cpp new file mode 100644 index 00000000000..7b8335ed909 --- /dev/null +++ b/src/KOKKOS/pair_granular_kokkos.cpp @@ -0,0 +1,406 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. +------------------------------------------------------------------------- */ + +#include "pair_granular_kokkos.h" + +#include "atom_kokkos.h" +#include "atom_masks.h" +#include "error.h" +#include "fix_neigh_history_kokkos.h" +#include "force.h" +#include "gran_sub_mod.h" +#include "granular_model.h" +#include "kokkos.h" +#include "modify.h" +#include "neigh_list_kokkos.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; +using namespace Granular_NS; + +namespace { +enum {NORMAL_HOOKE, NORMAL_HERTZ, NORMAL_HERTZ_MATERIAL}; +enum {DAMP_NONE, DAMP_VELOCITY, DAMP_MASS_VELOCITY, DAMP_VISCOELASTIC, DAMP_TSUJI, DAMP_COR}; +enum {TANG_NONE, TANG_LINEAR_NOHISTORY, TANG_LINEAR_HISTORY, TANG_MINDLIN, TANG_MINDLIN_FORCE}; +enum {ROLL_NONE, ROLL_SDS}; +enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; +enum {HEAT_NONE, HEAT_RADIUS, HEAT_AREA}; +} + +template +PairGranularKokkos::PairGranularKokkos(LAMMPS *lmp) : PairGranular(lmp) +{ + kokkosable = 1; + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + no_virial_fdotr_compute = 0; + datamask_read = X_MASK | V_MASK | OMEGA_MASK | F_MASK | TORQUE_MASK | TYPE_MASK | MASK_MASK | + RMASS_MASK | RADIUS_MASK; + datamask_modify = F_MASK | TORQUE_MASK; + fix_historyKK = nullptr; +} + +template +PairGranularKokkos::~PairGranularKokkos() +{ + if (copymode) return; +} + +template +const char *PairGranularKokkos::history_fix_style() const +{ + return execution_space == Device ? "NEIGH_HISTORY/KK/DEVICE" : "NEIGH_HISTORY/KK/HOST"; +} + +template +void PairGranularKokkos::validate_models() const +{ + const char *categories[] = {"normal", "damping", "tangential", "rolling", "twisting", "heat"}; + const char *supported[][6] = { + {"hooke", "hertz", "hertz/material", nullptr, nullptr, nullptr}, + {"none", "velocity", "mass_velocity", "viscoelastic", "tsuji", "coeff_restitution"}, + {"none", "linear_nohistory", "linear_history", "mindlin", "mindlin/force", nullptr}, + {"none", "sds", nullptr, nullptr, nullptr, nullptr}, + {"none", "sds", "marshall", nullptr, nullptr, nullptr}, + {"none", "radius", "area", nullptr, nullptr, nullptr}}; + + for (int n = 0; n < nmodels; ++n) { + auto *model = models_list[n]; + for (int category = 0; category < NSUBMODELS; ++category) { + const auto &name = model->sub_models[category]->name; + bool found = false; + for (int i = 0; supported[category][i]; ++i) + if (name == supported[category][i]) found = true; + if (!found) { + if ((name == "mindlin_rescale") || (name == "mindlin_rescale/force")) + error->all(FLERR, "Pair granular/kk does not support tangential model {} because it requires nondefault history transfer", name); + error->all(FLERR, "Pair granular/kk does not support {} model {}", categories[category], name); + } + } + if (model->nondefault_history_transfer) + error->all(FLERR, "Pair granular/kk does not support a model requiring nondefault history transfer"); + } +} + +template +void PairGranularKokkos::init_style() +{ + validate_models(); + PairGranular::init_style(); + + if (use_history) { + fix_historyKK = dynamic_cast *>(fix_history); + if (!fix_historyKK) error->all(FLERR, "Pair granular/kk requires NEIGH_HISTORY/KK"); + } + + const int neighflag = lmp->kokkos->neighflag; + if (neighflag == FULL) + error->all(FLERR, "Pair granular/kk requires a half neighbor list"); + auto *request = neighbor->find_request(this); + request->set_kokkos_host(std::is_same_v && + !std::is_same_v); + request->set_kokkos_device(std::is_same_v); +} + +template +double PairGranularKokkos::init_one(int i, int j) +{ + const double cutoff = PairGranular::init_one(i,j); + build_models(); + return cutoff; +} + +template +void PairGranularKokkos::build_models() +{ + const int ntypes = atom->ntypes; + k_models = Kokkos::DualView("granular:models", + ntypes + 1, ntypes + 1); + auto h_models = k_models.view_host(); + for (int i = 1; i <= ntypes; ++i) { + for (int j = 1; j <= ntypes; ++j) { + auto *gm = models_list[types_indices[i][j]]; + Model p{}; + auto *normal = gm->sub_models[NORMAL]; + auto *damping = gm->sub_models[DAMPING]; + auto *tangential = gm->sub_models[TANGENTIAL]; + auto *rolling = gm->sub_models[ROLLING]; + auto *twisting = gm->sub_models[TWISTING]; + auto *heat = gm->sub_models[HEAT]; + + if (normal->name == "hooke") { p.normal = NORMAL_HOOKE; p.normal_k = normal->coeffs[0]; } + else if (normal->name == "hertz") { p.normal = NORMAL_HERTZ; p.normal_k = normal->coeffs[0]; } + else { + p.normal = NORMAL_HERTZ_MATERIAL; + const auto e = normal->coeffs[0], nu = normal->coeffs[2]; + p.normal_k = static_cast(2.0 * e / (3.0 * (1.0 - nu * nu))); + } + const double normal_damp = normal->coeffs[1]; + if (damping->name == "none") { p.damping = DAMP_NONE; p.damping_coeff = 0; } + else if (damping->name == "velocity") { p.damping = DAMP_VELOCITY; p.damping_coeff = normal_damp; } + else if (damping->name == "mass_velocity") { p.damping = DAMP_MASS_VELOCITY; p.damping_coeff = normal_damp; } + else if (damping->name == "viscoelastic") { p.damping = DAMP_VISCOELASTIC; p.damping_coeff = normal_damp; } + else { + p.damping = damping->name == "tsuji" ? DAMP_TSUJI : DAMP_COR; + if (p.damping == DAMP_TSUJI) { + const double q = normal_damp; + p.damping_coeff = 1.2728 - 4.2783*q + 11.087*q*q - 22.348*q*q*q + + 27.467*q*q*q*q - 18.022*q*q*q*q*q + 4.8218*q*q*q*q*q*q; + } else { + const double l = std::log(normal_damp); + constexpr double pi = 3.14159265358979323846; + p.damping_coeff = normal->name == "hooke" ? + -2.0*l/std::sqrt(pi*pi+l*l) : + -1.8257418583505538*l/std::sqrt(pi*pi+l*l); + } + } + if (tangential->name == "none") p.tangential = TANG_NONE; + else if (tangential->name == "linear_nohistory") p.tangential = TANG_LINEAR_NOHISTORY; + else if (tangential->name == "linear_history") p.tangential = TANG_LINEAR_HISTORY; + else if (tangential->name == "mindlin") p.tangential = TANG_MINDLIN; + else p.tangential = TANG_MINDLIN_FORCE; + p.tangential_k = tangential->num_coeffs ? tangential->coeffs[0] : 0; + if ((p.tangential == TANG_MINDLIN) || (p.tangential == TANG_MINDLIN_FORCE)) { + if (p.tangential_k < 0) { + const double e = normal->coeffs[0], nu = normal->coeffs[2]; + p.tangential_k = 2.0 * e / ((2.0-nu)*(1.0+nu)); + } + } + p.tangential_xt = tangential->num_coeffs > 1 ? tangential->coeffs[1] : 0; + p.tangential_mu = tangential->num_coeffs > 1 ? + tangential->coeffs[tangential->num_coeffs - 1] : 0; + p.tangential_history = tangential->history_index; + + p.rolling = rolling->name == "sds" ? ROLL_SDS : ROLL_NONE; + p.rolling_k = rolling->num_coeffs ? rolling->coeffs[0] : 0; + p.rolling_damp = rolling->num_coeffs ? rolling->coeffs[1] : 0; + p.rolling_mu = rolling->num_coeffs ? rolling->coeffs[2] : 0; + p.rolling_history = rolling->history_index; + p.twisting = twisting->name == "sds" ? TWIST_SDS : + (twisting->name == "marshall" ? TWIST_MARSHALL : TWIST_NONE); + p.twisting_k = twisting->name == "marshall" ? p.tangential_k : (twisting->num_coeffs ? twisting->coeffs[0] : 0); + p.twisting_damp = twisting->name == "marshall" ? p.tangential_xt : (twisting->num_coeffs ? twisting->coeffs[1] : 0); + p.twisting_mu = twisting->name == "marshall" ? p.tangential_mu : (twisting->num_coeffs ? twisting->coeffs[2] : 0); + p.twisting_history = twisting->history_index; + p.heat = heat->name == "radius" ? HEAT_RADIUS : (heat->name == "area" ? HEAT_AREA : HEAT_NONE); + p.heat_coeff = heat->num_coeffs ? heat->coeffs[0] : 0; + p.limit_damping = gm->limit_damping; + p.synchronized_verlet = gm->synchronized_verlet; + h_models(i,j) = p; + } + } + k_models.modify_host(); + models = k_models.template view(); +} + +template +KOKKOS_INLINE_FUNCTION +void PairGranularKokkos::evaluate(const Model &p, Contact &c, Result &out, + const KK_FLOAT dt, const bool update) +{ + for (int k = 0; k < 3; ++k) out.force[k] = out.torquei[k] = out.torquej[k] = + out.fs[k] = out.fr[k] = 0; + out.heat = out.twist = 0; + const KK_FLOAT rsq = c.dx[0]*c.dx[0] + c.dx[1]*c.dx[1] + c.dx[2]*c.dx[2]; + const KK_FLOAT r = Kokkos::sqrt(rsq); + const KK_FLOAT rinv = 1.0/r; + const KK_FLOAT delta = c.radi + c.radj - r; + const KK_FLOAT reff = c.radi*c.radj/(c.radi+c.radj); + const KK_FLOAT a = Kokkos::sqrt(delta*reff); + KK_FLOAT n[3], vr[3], vt[3], vtr[3], wr[3]; + for (int k = 0; k < 3; ++k) { n[k] = c.dx[k]*rinv; vr[k] = c.vi[k]-c.vj[k]; } + const KK_FLOAT vnnr = vr[0]*n[0]+vr[1]*n[1]+vr[2]*n[2]; + for (int k = 0; k < 3; ++k) vt[k] = vr[k]-vnnr*n[k]; + for (int k = 0; k < 3; ++k) wr[k] = c.radi*c.wi[k]+c.radj*c.wj[k]; + vtr[0] = vt[0] - (wr[1]*n[2]-wr[2]*n[1]); + vtr[1] = vt[1] - (wr[2]*n[0]-wr[0]*n[2]); + vtr[2] = vt[2] - (wr[0]*n[1]-wr[1]*n[0]); + const KK_FLOAT vrel = Kokkos::sqrt(vtr[0]*vtr[0]+vtr[1]*vtr[1]+vtr[2]*vtr[2]); + const KK_FLOAT fnormal = p.normal == NORMAL_HOOKE ? p.normal_k*delta : p.normal_k*a*delta; + KK_FLOAT dpref = 0; + if (p.damping == DAMP_VELOCITY) dpref = p.damping_coeff; + else if (p.damping == DAMP_MASS_VELOCITY) dpref = p.damping_coeff*c.meff; + else if (p.damping == DAMP_VISCOELASTIC) dpref = p.damping_coeff*c.meff*a; + else if ((p.damping == DAMP_TSUJI) || (p.damping == DAMP_COR)) { + const KK_FLOAT q = delta > 0 ? c.meff*fnormal/delta : 0; + dpref = p.damping_coeff*Kokkos::sqrt(q > 0 ? q : 0); + } + KK_FLOAT fntot = fnormal-dpref*vnnr; + if (p.limit_damping && fntot < 0) fntot = 0; + const KK_FLOAT fncrit = Kokkos::fabs(fntot); + + if (p.tangential == TANG_LINEAR_NOHISTORY) { + const KK_FLOAT mag = p.tangential_xt*dpref*vrel; + const KK_FLOAT scale = vrel > 0 ? (mag < p.tangential_mu*fncrit ? mag : p.tangential_mu*fncrit)/vrel : 0; + for (int k = 0; k < 3; ++k) out.fs[k] = -scale*vtr[k]; + } else if (p.tangential != TANG_NONE) { + KK_FLOAT *h = c.history + p.tangential_history; + KK_FLOAT shear[3] = {h[0],h[1],h[2]}; + if (update) { + const KK_FLOAT dot = shear[0]*n[0]+shear[1]*n[1]+shear[2]*n[2]; + for (int k = 0; k < 3; ++k) shear[k] -= dot*n[k]; + const KK_FLOAT ks = p.tangential == TANG_LINEAR_HISTORY ? p.tangential_k : p.tangential_k*a; + for (int k = 0; k < 3; ++k) shear[k] += (p.tangential == TANG_MINDLIN_FORCE ? -ks*dt : dt)*vtr[k]; + } + const KK_FLOAT ks = p.tangential == TANG_LINEAR_HISTORY ? p.tangential_k : p.tangential_k*a; + for (int k = 0; k < 3; ++k) + out.fs[k] = p.tangential == TANG_MINDLIN_FORCE ? shear[k]-p.tangential_xt*dpref*vtr[k] : + -ks*shear[k]-p.tangential_xt*dpref*vtr[k]; + const KK_FLOAT mag = Kokkos::sqrt(out.fs[0]*out.fs[0]+out.fs[1]*out.fs[1]+out.fs[2]*out.fs[2]); + const KK_FLOAT lim = p.tangential_mu*fncrit; + if (mag > lim && mag > 0) for (int k = 0; k < 3; ++k) out.fs[k] *= lim/mag; + if (update) for (int k = 0; k < 3; ++k) h[k] = shear[k]; + } + for (int k = 0; k < 3; ++k) out.force[k] = fntot*n[k]+out.fs[k]; + const KK_FLOAT di = c.radi-0.5*delta, dj = c.radj-0.5*delta; + const KK_FLOAT cross[3] = {n[1]*out.fs[2]-n[2]*out.fs[1], n[2]*out.fs[0]-n[0]*out.fs[2], + n[0]*out.fs[1]-n[1]*out.fs[0]}; + for (int k = 0; k < 3; ++k) { out.torquei[k] = -di*cross[k]; out.torquej[k] = -dj*cross[k]; } + + if (p.rolling == ROLL_SDS) { + KK_FLOAT *h = c.history+p.rolling_history, rel[3], vrl[3], fr[3]; + for (int k = 0; k < 3; ++k) rel[k] = c.wi[k]-c.wj[k]; + vrl[0] = reff*(rel[1]*n[2]-rel[2]*n[1]); vrl[1] = reff*(rel[2]*n[0]-rel[0]*n[2]); + vrl[2] = reff*(rel[0]*n[1]-rel[1]*n[0]); + if (update) for (int k = 0; k < 3; ++k) h[k] += dt*vrl[k]; + for (int k = 0; k < 3; ++k) fr[k] = -p.rolling_k*h[k]-p.rolling_damp*vrl[k]; + const KK_FLOAT mag = Kokkos::sqrt(fr[0]*fr[0]+fr[1]*fr[1]+fr[2]*fr[2]), lim = p.rolling_mu*fncrit; + if (mag > lim && mag > 0) for (int k = 0; k < 3; ++k) fr[k] *= lim/mag; + const KK_FLOAT tr[3] = {reff*(n[1]*fr[2]-n[2]*fr[1]), reff*(n[2]*fr[0]-n[0]*fr[2]), + reff*(n[0]*fr[1]-n[1]*fr[0])}; + for (int k = 0; k < 3; ++k) { out.fr[k] = fr[k]; out.torquei[k] += tr[k]; out.torquej[k] -= tr[k]; } + } + if (p.twisting != TWIST_NONE) { + KK_FLOAT *h = c.history+p.twisting_history; + const KK_FLOAT twist = (c.wi[0]-c.wj[0])*n[0]+(c.wi[1]-c.wj[1])*n[1]+(c.wi[2]-c.wj[2])*n[2]; + const KK_FLOAT k = p.twisting == TWIST_MARSHALL ? 0.5*p.twisting_k*a*a : p.twisting_k; + const KK_FLOAT damp = p.twisting == TWIST_MARSHALL ? 0.5*p.twisting_damp*dpref*a*a : p.twisting_damp; + const KK_FLOAT lim = (p.twisting == TWIST_MARSHALL ? 2.0*a/3.0*p.twisting_mu : p.twisting_mu)*fncrit; + if (update) h[0] += twist*dt; + KK_FLOAT mt = -k*h[0]-damp*twist; + if (Kokkos::fabs(mt) > lim) mt = mt > 0 ? lim : -lim; + out.twist = mt; + for (int q = 0; q < 3; ++q) { out.torquei[q] += mt*n[q]; out.torquej[q] -= mt*n[q]; } + } + if (p.heat == HEAT_RADIUS) out.heat = 2*p.heat_coeff*a*(c.tj-c.ti); + else if (p.heat == HEAT_AREA) out.heat = p.heat_coeff*3.141592653589793*a*a*(c.tj-c.ti); +} + +template +KOKKOS_INLINE_FUNCTION +void PairGranularKokkos::operator()(TagPairGranularKokkosCompute, const int ii) const +{ + const int i = d_ilist[ii]; + const int jnum = d_numneigh[i]; + for (int jj = 0; jj < jnum; ++jj) { + int j = d_neighbors(i,jj); + const KK_FLOAT factor = special_lj[sbmask(j)]; + j &= NEIGHMASK; + if (factor == 0) continue; + Contact c{}; + for (int k = 0; k < 3; ++k) { c.dx[k] = x(i,k)-x(j,k); c.vi[k] = v(i,k); c.vj[k] = v(j,k); + c.wi[k] = omega(i,k); c.wj[k] = omega(j,k); } + c.radi = radius[i]; c.radj = radius[j]; + const KK_FLOAT rsq = c.dx[0]*c.dx[0]+c.dx[1]*c.dx[1]+c.dx[2]*c.dx[2]; + if (rsq >= (c.radi+c.radj)*(c.radi+c.radj)) { + if (use_history) { d_firsttouch(i,jj) = 0; for (int k = 0; k < size_history; ++k) d_firsthistory(i,size_history*jj+k) = 0; } + continue; + } + const KK_FLOAT mi = rmass[i], mj = rmass[j]; + c.meff = mi*mj/(mi+mj); + if (mask[i] & freeze_group_bit) c.meff = mj; + if (mask[j] & freeze_group_bit) c.meff = mi; + c.ti = heat_flag ? d_temperature[i] : 0; c.tj = heat_flag ? d_temperature[j] : 0; + c.history = use_history ? &d_firsthistory(i,size_history*jj) : nullptr; + if (use_history) d_firsttouch(i,jj) = 1; + Result out{}; + evaluate(models(type[i],type[j]),c,out,dt_kk,history_update); + for (int k = 0; k < 3; ++k) { + Kokkos::atomic_add(&f(i,k),factor*out.force[k]); Kokkos::atomic_add(&torque(i,k),factor*out.torquei[k]); + if (newton_pair || j < nlocal) { Kokkos::atomic_add(&f(j,k),-factor*out.force[k]); + Kokkos::atomic_add(&torque(j,k),factor*out.torquej[k]); } + } + if (heat_flag) { Kokkos::atomic_add(&d_heatflow[i],out.heat); + if (newton_pair || j < nlocal) Kokkos::atomic_add(&d_heatflow[j],-out.heat); } + } +} + +template +void PairGranularKokkos::compute(int eflag_in, int vflag_in) +{ + ev_init(eflag_in,vflag_in,0); + atomKK->sync(execution_space,datamask_read); + atomKK->modified(execution_space,datamask_modify); + k_models.template sync(); + x = atomKK->k_x.view(); v = atomKK->k_v.view(); + omega = atomKK->k_omega.view(); f = atomKK->k_f.view(); + torque = atomKK->k_torque.view(); type = atomKK->k_type.view(); + mask = atomKK->k_mask.view(); rmass = atomKK->k_rmass.view(); + radius = atomKK->k_radius.view(); + nlocal = atom->nlocal; newton_pair = force->newton_pair; history_update = update->setupflag == 0; dt_kk = update->dt; + for (int k = 0; k < 4; ++k) special_lj[k] = force->special_lj[k]; + auto *k_list = static_cast *>(list); + d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; d_numneigh = k_list->d_numneigh; + if (use_history) { + fix_historyKK->k_firstflag.template sync(); fix_historyKK->k_firstvalue.template sync(); + d_firsttouch = fix_historyKK->k_firstflag.template view(); + d_firsthistory = fix_historyKK->k_firstvalue.template view(); + } + if (heat_flag) { + const int nall = atom->nlocal + atom->nghost; + k_temperature = Kokkos::DualView("granular:temperature",nall); + k_heatflow = Kokkos::DualView("granular:heatflow",nall); + auto ht = k_temperature.view_host(), hh = k_heatflow.view_host(); + for (int i = 0; i < nall; ++i) { ht(i) = atom->temperature[i]; hh(i) = 0; } + k_temperature.modify_host(); k_heatflow.modify_host(); k_temperature.template sync(); + k_heatflow.template sync(); d_temperature = k_temperature.template view(); d_heatflow = k_heatflow.template view(); + } + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,list->inum),*this); + copymode = 0; + if (use_history) { fix_historyKK->k_firstflag.template modify(); fix_historyKK->k_firstvalue.template modify(); } + if (heat_flag) { + k_heatflow.template modify(); k_heatflow.sync_host(); + auto hh = k_heatflow.view_host(); + for (int i = 0; i < nlocal; ++i) atom->heatflow[i] += hh(i); + } + if (vflag_fdotr) pair_virial_fdotr_compute(this); +} + +template +double PairGranularKokkos::single(int i, int j, int itype, int jtype, + double, double, double, double &fforce) +{ + Model p = k_models.view_host()(itype,jtype); + Contact c{}; + for (int k = 0; k < 3; ++k) { c.dx[k] = atom->x[i][k]-atom->x[j][k]; c.vi[k] = atom->v[i][k]; + c.vj[k] = atom->v[j][k]; c.wi[k] = atom->omega[i][k]; c.wj[k] = atom->omega[j][k]; } + c.radi = atom->radius[i]; c.radj = atom->radius[j]; c.meff = atom->rmass[i]*atom->rmass[j]/(atom->rmass[i]+atom->rmass[j]); + KK_FLOAT history[9] = {}; c.history = history; + Result out{}; evaluate(p,c,out,update->dt,false); + const double r = std::sqrt(c.dx[0]*c.dx[0]+c.dx[1]*c.dx[1]+c.dx[2]*c.dx[2]); + fforce = (out.force[0]*c.dx[0]+out.force[1]*c.dx[1]+out.force[2]*c.dx[2])/(r*r); + for (int k = 0; k < single_extra; ++k) svector[k] = 0; + for (int k = 0; k < 3; ++k) { svector[k] = out.fs[k]; svector[4+k] = out.fr[k]; svector[9+k] = c.dx[k]; } + svector[3] = std::sqrt(out.fs[0]*out.fs[0]+out.fs[1]*out.fs[1]+out.fs[2]*out.fs[2]); + svector[7] = std::sqrt(out.fr[0]*out.fr[0]+out.fr[1]*out.fr[1]+out.fr[2]*out.fr[2]); svector[8] = out.twist; + return 0.0; +} + +template class PairGranularKokkos; +#ifdef LMP_KOKKOS_GPU +template class PairGranularKokkos; +#endif diff --git a/src/KOKKOS/pair_granular_kokkos.h b/src/KOKKOS/pair_granular_kokkos.h new file mode 100644 index 00000000000..924b081fb1e --- /dev/null +++ b/src/KOKKOS/pair_granular_kokkos.h @@ -0,0 +1,117 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(granular/kk,PairGranularKokkos); +PairStyle(granular/kk/device,PairGranularKokkos); +PairStyle(granular/kk/host,PairGranularKokkos); +// clang-format on +#else + +#ifndef LMP_PAIR_GRANULAR_KOKKOS_H +#define LMP_PAIR_GRANULAR_KOKKOS_H + +#include "kokkos_type.h" +#include "pair_granular.h" +#include "pair_kokkos.h" + +namespace LAMMPS_NS { + +template class FixNeighHistoryKokkos; +struct TagPairGranularKokkosCompute {}; + +template +class PairGranularKokkos : public PairGranular { + public: + typedef DeviceType device_type; + typedef ArrayTypes AT; + + PairGranularKokkos(class LAMMPS *); + ~PairGranularKokkos() override; + void compute(int, int) override; + void init_style() override; + double init_one(int, int) override; + double single(int, int, int, int, double, double, double, double &) override; + + struct Model { + int normal, damping, tangential, rolling, twisting, heat; + int tangential_history, rolling_history, twisting_history; + int limit_damping, synchronized_verlet; + KK_FLOAT normal_k, damping_coeff; + KK_FLOAT tangential_k, tangential_xt, tangential_mu; + KK_FLOAT rolling_k, rolling_damp, rolling_mu; + KK_FLOAT twisting_k, twisting_damp, twisting_mu; + KK_FLOAT heat_coeff; + }; + + protected: + typename AT::t_kkfloat_1d_3_lr_randomread x; + typename AT::t_kkfloat_1d_3_randomread v; + typename AT::t_kkfloat_1d_3_randomread omega; + typename AT::t_kkacc_1d_3 f; + typename AT::t_kkacc_1d_3 torque; + typename AT::t_int_1d_randomread type; + typename AT::t_int_1d_randomread mask; + typename AT::t_kkfloat_1d_randomread rmass; + typename AT::t_kkfloat_1d_randomread radius; + typename AT::t_neighbors_2d d_neighbors; + typename AT::t_int_1d_randomread d_ilist; + typename AT::t_int_1d_randomread d_numneigh; + typename AT::t_int_2d d_firsttouch; + typename AT::t_kkfloat_2d d_firsthistory; + + Kokkos::DualView k_models; + typename Kokkos::DualView::t_dev_const models; + Kokkos::DualView k_temperature, k_heatflow; + typename AT::t_kkfloat_1d d_temperature, d_heatflow; + + FixNeighHistoryKokkos *fix_historyKK; + int nlocal, newton_pair, history_update; + KK_FLOAT dt_kk; + KK_FLOAT special_lj[4]; + + const char *history_fix_style() const override; + void build_models(); + void validate_models() const; + + struct Contact { + KK_FLOAT dx[3], vi[3], vj[3], wi[3], wj[3]; + KK_FLOAT radi, radj, meff, ti, tj; + KK_FLOAT *history; + }; + struct Result { + KK_FLOAT force[3], torquei[3], torquej[3], heat; + KK_FLOAT fs[3], fr[3], twist; + }; + +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static void evaluate(const Model &, Contact &, Result &, const KK_FLOAT, const bool); + + public: +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(TagPairGranularKokkosCompute, const int) const; + +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static int sbmask(const int j) { return j >> SBBITS & 3; } + + friend void pair_virial_fdotr_compute(PairGranularKokkos *); +}; + +} // namespace LAMMPS_NS + +#endif +#endif