diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.cpp b/src/KOKKOS/fix_wall_lj93_kokkos.cpp index 93bcc9abc9..fe75c92a8a 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj93_kokkos.cpp @@ -12,13 +12,26 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Mitch Murphy (alphataubio@gmail.com) +------------------------------------------------------------------------- */ + #include "fix_wall_lj93_kokkos.h" +#include "atom_masks.h" #include "atom_kokkos.h" #include "error.h" -#include "atom_masks.h" +#include "input.h" +#include "math_special.h" +#include "memory_kokkos.h" +#include "modify_kokkos.h" +#include "variable.h" +#include "update.h" + +#include using namespace LAMMPS_NS; +using MathSpecial::powint; /* ---------------------------------------------------------------------- */ @@ -29,11 +42,141 @@ FixWallLJ93Kokkos::FixWallLJ93Kokkos(LAMMPS *lmp, int narg, char **a kokkosable = 1; atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; - datamask_read = EMPTY_MASK; - datamask_modify = EMPTY_MASK; + datamask_read = X_MASK | V_MASK | MASK_MASK; + datamask_modify = F_MASK; virial_global_flag = virial_peratom_flag = 0; + + memoryKK->create_kokkos(k_epsilon,6,"wall_lj93:epsilon"); + memoryKK->create_kokkos(k_sigma,6,"wall_lj93:sigma"); + memoryKK->create_kokkos(k_cutoff,6,"wall_lj93:cutoff"); + + d_epsilon = k_epsilon.template view(); + d_sigma = k_sigma.template view(); + d_cutoff = k_cutoff.template view(); + + for( int i=0 ; i<6 ; i++ ) { + k_epsilon.h_view(i) = epsilon[i]; + k_sigma.h_view(i) = sigma[i]; + k_cutoff.h_view(i) = cutoff[i]; + } + + k_epsilon.template modify(); + k_sigma.template modify(); + k_cutoff.template modify(); + + k_epsilon.template sync(); + k_sigma.template sync(); + k_cutoff.template sync(); + + memoryKK->create_kokkos(d_coeff1,6,"wall_lj93:coeff1"); + memoryKK->create_kokkos(d_coeff2,6,"wall_lj93:coeff2"); + memoryKK->create_kokkos(d_coeff3,6,"wall_lj93:coeff3"); + memoryKK->create_kokkos(d_coeff4,6,"wall_lj93:coeff4"); + memoryKK->create_kokkos(d_offset,6,"wall_lj93:offset"); + + memoryKK->create_kokkos(k_ewall,ewall,7,"wall_lj93:ewall"); + d_ewall = k_ewall.template view(); + } +template +FixWallLJ93Kokkos::~FixWallLJ93Kokkos() +{ + if (copymode) return; + + memoryKK->destroy_kokkos(k_epsilon); + memoryKK->destroy_kokkos(k_sigma); + memoryKK->destroy_kokkos(k_cutoff); + + memoryKK->destroy_kokkos(d_coeff1); + memoryKK->destroy_kokkos(d_coeff2); + memoryKK->destroy_kokkos(d_coeff3); + memoryKK->destroy_kokkos(d_coeff4); + memoryKK->destroy_kokkos(d_offset); + + //std::cerr << "ok 2\n"; + +} + +/* ---------------------------------------------------------------------- */ + +template +void FixWallLJ93Kokkos::precompute(int m) +{ + d_coeff1[m] = 6.0 / 5.0 * d_epsilon[m] * powint(d_sigma[m], 9); + d_coeff2[m] = 3.0 * d_epsilon[m] * powint(d_sigma[m], 3); + d_coeff3[m] = 2.0 / 15.0 * d_epsilon[m] * powint(d_sigma[m], 9); + d_coeff4[m] = d_epsilon[m] * powint(d_sigma[m], 3); + + double rinv = 1.0 / d_cutoff[m]; + double r2inv = rinv * rinv; + double r4inv = r2inv * r2inv; + d_offset[m] = d_coeff3[m] * r4inv * r4inv * rinv - d_coeff4[m] * r2inv * rinv; +} + + +/* ---------------------------------------------------------------------- */ + +template +void FixWallLJ93Kokkos::post_force(int vflag) +{ + + //std::cerr << "post_force DeviceType=" << DeviceType << "\n"; + + atomKK->sync(execution_space,datamask_read); + atomKK->modified(execution_space,datamask_modify); + + // virial setup + + v_init(vflag); + + // energy intialize. + // eflag is used to track whether wall energies have been communicated. + + eflag = 0; + //for (int m = 0; m <= nwall; m++) d_ewall(m) = 0.0; + for (int m = 0; m <= nwall; m++) k_ewall.d_view(m) = 0.0; + + // coord = current position of wall + // evaluate variables if necessary, wrap with clear/add + // for epsilon/sigma variables need to re-invoke precompute() + + if (varflag) modify->clearstep_compute(); + + double coord; + for (int m = 0; m < nwall; m++) { + if (xstyle[m] == VARIABLE) { + coord = input->variable->compute_equal(xindex[m]); + if (wallwhich[m] < YLO) + coord *= xscale; + else if (wallwhich[m] < ZLO) + coord *= yscale; + else + coord *= zscale; + } else + coord = coord0[m]; + if (wstyle[m] == VARIABLE) { + if (estyle[m] == VARIABLE) { + d_epsilon[m] = input->variable->compute_equal(eindex[m]); + if (d_epsilon[m] < 0.0) error->all(FLERR, "Variable evaluation in fix wall gave bad value"); + } + if (sstyle[m] == VARIABLE) { + d_sigma[m] = input->variable->compute_equal(sindex[m]); + if (d_sigma[m] < 0.0) error->all(FLERR, "Variable evaluation in fix wall gave bad value"); + } + precompute(m); + } + + wall_particle(m, wallwhich[m], coord); + } + + k_ewall.template modify(); + k_ewall.template sync(); + + if (varflag) modify->addstep_compute(update->ntimestep + 1); +} + + /* ---------------------------------------------------------------------- interaction of all particles in group with a wall m = index of wall coeffs @@ -42,17 +185,16 @@ FixWallLJ93Kokkos::FixWallLJ93Kokkos(LAMMPS *lmp, int narg, char **a ------------------------------------------------------------------------- */ template +KOKKOS_INLINE_FUNCTION void FixWallLJ93Kokkos::wall_particle(int m_in, int which, double coord_in) { m = m_in; coord = coord_in; - - atomKK->sync(execution_space, X_MASK|F_MASK|MASK_MASK); - x = atomKK->k_x.view(); - f = atomKK->k_f.view(); - mask = atomKK->k_mask.view(); - - int nlocal = atom->nlocal; + d_x = atomKK->k_x.template view(); + d_f = atomKK->k_f.template view(); + d_mask = atomKK->k_mask.template view(); + int nlocal = atomKK->nlocal; + double tmp[7]; dim = which / 2; side = which % 2; @@ -60,34 +202,68 @@ void FixWallLJ93Kokkos::wall_particle(int m_in, int which, double co copymode = 1; FixWallLJ93KokkosFunctor wp_functor(this); - Kokkos::parallel_reduce(nlocal,wp_functor,ewall); + Kokkos::parallel_reduce(nlocal,wp_functor,tmp); copymode = 0; - atomKK->modified(execution_space, F_MASK); + //std::cerr << fmt::format("tmp[0]={} tmp[{}]={} \n",tmp[0],m+1,tmp[m+1]); + + + Kokkos::atomic_add(&(d_ewall[0]),tmp[0]); + Kokkos::atomic_add(&(d_ewall[m+1]),tmp[m+1]); + + //std::cerr << fmt::format("k_ewall.d_view[0]={} k_ewall.d_view[{}]={} \n",k_ewall.d_view[0],m+1,k_ewall.d_view[m+1]); + + } template KOKKOS_INLINE_FUNCTION void FixWallLJ93Kokkos::wall_particle_item(int i, value_type ewall) const { - if (mask(i) & groupbit) { + if (d_mask(i) & groupbit) { double delta; - if (side < 0) delta = x(i,dim) - coord; - else delta = coord - x(i,dim); - if (delta >= cutoff[m]) return; + if (side < 0) delta = d_x(i,dim) - coord; + else delta = coord - d_x(i,dim); + if (delta >= d_cutoff(m)) return; if (delta <= 0.0) Kokkos::abort("Particle on or inside fix wall surface"); double rinv = 1.0/delta; double r2inv = rinv*rinv; double r4inv = r2inv*r2inv; double r10inv = r4inv*r4inv*r2inv; - double fwall = side * (coeff1[m]*r10inv - coeff2[m]*r4inv); - f(i,dim) -= fwall; - ewall[0] += coeff3[m]*r4inv*r4inv*rinv - - coeff4[m]*r2inv*rinv - offset[m]; + double fwall = side * (d_coeff1(m)*r10inv - d_coeff2(m)*r4inv); + d_f(i,dim) -= fwall; + ewall[0] += d_coeff3(m)*r4inv*r4inv*rinv - d_coeff4(m)*r2inv*rinv - d_offset(m); ewall[m+1] += fwall; } + } +/* ---------------------------------------------------------------------- + energy of wall interaction +------------------------------------------------------------------------- */ + +/* +template +double FixWallLJ93Kokkos::compute_scalar() +{ + // only sum across procs one time + + //std::cerr << fmt::format("k_ewall[0] = {} d_ewall[0] = {}\n", k_ewall.h_view[0], k_ewall.d_view[0] ); + + k_ewall.template sync(); + + if (eflag == 0) { + MPI_Allreduce(k_ewall.h_view.data(), ewall_all, 7, MPI_DOUBLE, MPI_SUM, world); + eflag = 1; + } + + std::cerr << fmt::format("compute_scalar() = {}\n", ewall_all[0]); + return ewall_all[0]; + +} +*/ + + namespace LAMMPS_NS { template class FixWallLJ93Kokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.h b/src/KOKKOS/fix_wall_lj93_kokkos.h index 720e586f5d..26b50052fd 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.h +++ b/src/KOKKOS/fix_wall_lj93_kokkos.h @@ -36,8 +36,11 @@ class FixWallLJ93Kokkos : public FixWallLJ93 { typedef double value_type[]; FixWallLJ93Kokkos(class LAMMPS *, int, char **); + ~FixWallLJ93Kokkos() override; + void precompute(int) override; + void post_force(int) override; void wall_particle(int, int, double) override; - + int m; KOKKOS_INLINE_FUNCTION @@ -47,9 +50,18 @@ class FixWallLJ93Kokkos : public FixWallLJ93 { int dim,side; double coord; - typename AT::t_x_array x; - typename AT::t_f_array f; - typename AT::t_int_1d mask; + typename AT::t_x_array d_x; + typename AT::t_f_array d_f; + typename AT::t_int_1d d_mask; + + typename AT::tdual_ffloat_1d k_epsilon,k_sigma,k_cutoff; + typename AT::t_ffloat_1d d_epsilon,d_sigma,d_cutoff; + + typename AT::t_ffloat_1d d_coeff1,d_coeff2,d_coeff3,d_coeff4,d_offset; + + typename AT::tdual_ffloat_1d k_ewall; + typename AT::t_ffloat_1d d_ewall; + }; template @@ -60,13 +72,16 @@ struct FixWallLJ93KokkosFunctor { FixWallLJ93Kokkos c; FixWallLJ93KokkosFunctor(FixWallLJ93Kokkos* c_ptr): - value_count(c_ptr->m+1), c(*c_ptr) {} + //value_count(c_ptr->m+1), c(*c_ptr) {} + value_count(7), c(*c_ptr) {} + KOKKOS_INLINE_FUNCTION void operator()(const int i, value_type ewall) const { c.wall_particle_item(i,ewall); } }; + } #endif diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 50289d0f69..b8fe84ea91 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -27,8 +27,6 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum { XLO = 0, XHI = 1, YLO = 2, YHI = 3, ZLO = 4, ZHI = 5 }; - static const char *wallpos[] = {"xlo", "xhi", "ylo", "yhi", "zlo", "zhi"}; /* ---------------------------------------------------------------------- */ @@ -44,6 +42,7 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall virial_global_flag = virial_peratom_flag = 1; respa_level_support = 1; ilevel_respa = 0; + ewall = new double[7]; // parse args diff --git a/src/fix_wall.h b/src/fix_wall.h index 81abfab8ea..85c37f8f17 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -20,6 +20,7 @@ namespace LAMMPS_NS { class FixWall : public Fix { public: + enum { XLO = 0, XHI = 1, YLO = 2, YHI = 3, ZLO = 4, ZHI = 5 }; int nwall; int wallwhich[6]; double coord0[6]; @@ -47,7 +48,7 @@ class FixWall : public Fix { protected: double epsilon[6], sigma[6], alpha[6], cutoff[6]; - double ewall[7], ewall_all[7]; + double *ewall, ewall_all[7]; // need ewall double*, not double[] for kokkos dual view double xscale, yscale, zscale; int estyle[6], sstyle[6], astyle[6], wstyle[6]; int eindex[6], sindex[6]; diff --git a/unittest/force-styles/tests/fix-timestep-wall_lj93_const.yaml b/unittest/force-styles/tests/fix-timestep-wall_lj93_const.yaml index c574eddf76..1756751e5e 100644 --- a/unittest/force-styles/tests/fix-timestep-wall_lj93_const.yaml +++ b/unittest/force-styles/tests/fix-timestep-wall_lj93_const.yaml @@ -11,7 +11,7 @@ pre_commands: ! | post_commands: ! | fix move all nve fix test solute wall/lj93 ylo EDGE 100.0 2.0 5.0 yhi EDGE 100.0 2.0 5.0 - fix_modify test virial yes +# fix_modify test virial yes input_file: in.fourmol natoms: 29 run_stress: ! |2-