diff --git a/doc/src/compute_stress_atom.rst b/doc/src/compute_stress_atom.rst index 6c4e0b690c7..15e5fe4539a 100644 --- a/doc/src/compute_stress_atom.rst +++ b/doc/src/compute_stress_atom.rst @@ -276,6 +276,8 @@ i.e. the kspace_style command in LAMMPS. It also does not implement the following fixes which add rigid-body constraints: :doc:`fix rigid/* ` and the OpenMP accelerated version of :doc:`fix rigid/small `, while all other :doc:`fix rigid/*/small ` are implemented. +The KOKKOS version of :doc:`fix shake ` is not implemented +either, while the non-accelerated version is. LAMMPS will generate an error if one of these options is included in your model. Extension of centroid stress calculations to these force diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index 268135d945a..0ef293a8ee8 100644 --- a/doc/src/fix_shake.rst +++ b/doc/src/fix_shake.rst @@ -304,6 +304,14 @@ can make minimization very inefficient and also cause stability problems with some minimization algorithms. Sometimes those can be avoided by reducing the :doc:`timestep `. +.. versionchanged:: TBD + +The *shake/kk* style does not compute the centroid virial of the +constraint forces, so it cannot be used with :doc:`compute +centroid/stress/atom `. Use the non-accelerated +*shake* style for that. Previously this combination was accepted but +silently contributed zero. + Related commands """""""""""""""" diff --git a/doc/src/pair_pace.rst b/doc/src/pair_pace.rst index 04060afbb2d..90aa8c39461 100644 --- a/doc/src/pair_pace.rst +++ b/doc/src/pair_pace.rst @@ -24,10 +24,14 @@ Syntax .. parsed-literal:: - keyword = *product* or *recursive* or *chunksize* + keyword = *product* or *recursive* or *chunksize* or *neigh* *product* = use product algorithm for basis functions *recursive* = use recursive algorithm for basis functions *chunksize* value = number of atoms in each pass + *neigh* value = *auto* or *shared* or *global* + *auto* = automatically select team scratch memory level for the neighbor list build (default) + *shared* = force on-chip (level 0) shared memory scratch + *global* = force global (level 1) memory scratch .. code-block:: LAMMPS @@ -90,6 +94,30 @@ For example if there are 8192 atoms in the simulation and the *chunksize* is set to 4096, the ACE calculation will be broken up into two passes (running on a single GPU). +.. versionadded:: TBD + +The keyword *neigh* is only recognized by the KOKKOS versions of the pair +styles (*pace/kk* and *pace/extrapolation/kk*); the non-accelerated styles +stop with an "unknown keyword" error if it is given. It only has an effect +on GPU backends. It controls which level of Kokkos team scratch memory is +used to build the short neighbor list. Level 0 is fast on-chip shared +memory, but it is a limited resource that can be exceeded when atoms +have many neighbors and/or when there are many atomic species, which +would otherwise abort the run with an error such as "Requested too much +scratch memory on level 0". Level 1 is (much larger) global memory, +which avoids the limit at the cost of slower access. + +With the default value *auto*, the pair style queries the amount of +shared memory available on the device (rather than assuming a fixed +value such as 48 KiB, so that larger limits available in newer versions +of the Kokkos library are used automatically) and transparently falls +back to level 1 when the request does not fit into level 0, printing a +warning the first time this happens. The value *shared* forces the use +of level 0 (on-chip) scratch memory, and *global* forces the use of +level 1 (global) scratch memory; the latter can be used to silence the +fallback warning or to force global memory when the automatic heuristic +is too conservative. + Extrapolation grade """"""""""""""""""" diff --git a/doc/src/run_style.rst b/doc/src/run_style.rst index af43a9b244d..f95e8928901 100644 --- a/doc/src/run_style.rst +++ b/doc/src/run_style.rst @@ -375,6 +375,15 @@ REPLICA package. Correspondingly the *respa/omp* style is available only if the OPENMP package was included. See the :doc:`Build package ` page for more info. +.. versionchanged:: TBD + +The *respa* style cannot be used together with the KOKKOS package. It +keeps a copy of the forces of each level and clears and sums those +through the plain LAMMPS arrays, and it calls the force computations +without the transfers between the host and the device that the *verlet* +style performs for KOKKOS. On a GPU this gives wrong forces with no +indication that anything is amiss, so the combination is refused. + Run style *verlet/split* is not compatible with kspace styles from the INTEL package and it is not compatible with any tip4p, dipole, or spin kspace styles. diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp index f6910826869..f34ce8d8556 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -19,6 +19,7 @@ #include "angle.h" #include "atom.h" +#include "atom_masks.h" #include "atom_vec.h" #include "bond.h" #include "comm.h" @@ -481,6 +482,8 @@ void ComputeBornMatrix::compute_numdiff() // store copy of current forces for owned and ghost atoms + atom->sync_host_arrays(X_MASK | F_MASK); + double **x = atom->x; double **f = atom->f; @@ -534,8 +537,12 @@ void ComputeBornMatrix::compute_numdiff() // restore original forces for owned and ghost atoms + atom->sync_host_arrays(F_MASK); + for (int i = 0; i < nall; i++) for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k]; + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- @@ -544,6 +551,12 @@ void ComputeBornMatrix::compute_numdiff() void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) { + // the strain goes into the plain coordinate array, and the energy evaluation + // that follows works from the KOKKOS copies, so bring the host side up to + // date first and hand the write over afterwards + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; // NOTE: virial_addon() expressions predicated on @@ -564,6 +577,8 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) x[i][k] = temp_x[i][k] + 0.5 * numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); x[i][l] = temp_x[i][l] + 0.5 * numdelta * magnitude * (temp_x[i][k] - fixedpoint[k]); } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -577,6 +592,9 @@ void ComputeBornMatrix::restore_atoms(int nall, int idir) int k = dirlist[idir][0]; int l = dirlist[idir][1]; + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; if (l == k) for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k]; @@ -585,6 +603,8 @@ void ComputeBornMatrix::restore_atoms(int nall, int idir) x[i][l] = temp_x[i][l]; x[i][k] = temp_x[i][k]; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -674,9 +694,16 @@ void ComputeBornMatrix::virial_addon() void ComputeBornMatrix::force_clear(int nall) { + // the forces are cleared through the plain array and accumulated again by + // the force computations, which work from the KOKKOS copies + + atom->sync_host_arrays(F_MASK); + double **forces = atom->f; size_t nbytes = 3 * sizeof(double) * nall; if (nbytes) memset(&forces[0][0], 0, nbytes); + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-FIX/fix_numdiff.cpp b/src/EXTRA-FIX/fix_numdiff.cpp index fbebe159e61..acd7e222e24 100644 --- a/src/EXTRA-FIX/fix_numdiff.cpp +++ b/src/EXTRA-FIX/fix_numdiff.cpp @@ -19,6 +19,7 @@ #include "angle.h" #include "atom.h" +#include "atom_masks.h" #include "bond.h" #include "compute.h" #include "dihedral.h" @@ -191,6 +192,8 @@ void FixNumDiff::calculate_forces() // store copy of current forces for owned and ghost atoms + atom->sync_host_arrays(X_MASK | F_MASK); + double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; @@ -245,8 +248,12 @@ void FixNumDiff::calculate_forces() // restore original forces for owned and ghost atoms + atom->sync_host_arrays(F_MASK); + for (i = 0; i < nall; i++) for (j = 0; j < 3; j++) { f[i][j] = temp_f[i][j]; } + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- @@ -257,6 +264,12 @@ void FixNumDiff::displace_atoms(int ilocal, int idim, int magnitude) { if (ilocal < 0) return; + // the displacement goes into the plain coordinate array, and the energy + // evaluation that follows works from the KOKKOS copies, so bring the host + // side up to date first and hand the write over afterwards + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int *sametag = atom->sametag; int j = ilocal; @@ -266,6 +279,8 @@ void FixNumDiff::displace_atoms(int ilocal, int idim, int magnitude) j = sametag[j]; x[j][idim] += delta * magnitude; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -276,6 +291,8 @@ void FixNumDiff::restore_atoms(int ilocal, int idim) { if (ilocal < 0) return; + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int *sametag = atom->sametag; int j = ilocal; @@ -285,6 +302,8 @@ void FixNumDiff::restore_atoms(int ilocal, int idim) j = sametag[j]; x[j][idim] = temp_x[j][idim]; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -294,7 +313,12 @@ void FixNumDiff::restore_atoms(int ilocal, int idim) double FixNumDiff::update_energy() { + // the forces are cleared through the plain array and accumulated again by + // the force computations, which work from the KOKKOS copies + + atom->sync_host_arrays(F_MASK); force_clear(atom->f); + atom->modified_host_arrays(F_MASK); // flag that we only need to compute the global energy int eflag = ENERGY_GLOBAL | ENERGY_ONLY; diff --git a/src/EXTRA-FIX/fix_numdiff_virial.cpp b/src/EXTRA-FIX/fix_numdiff_virial.cpp index beff4b616d2..e6b8ba14ef4 100644 --- a/src/EXTRA-FIX/fix_numdiff_virial.cpp +++ b/src/EXTRA-FIX/fix_numdiff_virial.cpp @@ -19,6 +19,7 @@ #include "angle.h" #include "atom.h" +#include "atom_masks.h" #include "bond.h" #include "compute.h" #include "dihedral.h" @@ -201,6 +202,8 @@ void FixNumDiffVirial::calculate_virial() // store copy of current forces for owned and ghost atoms + atom->sync_host_arrays(X_MASK | F_MASK); + double **x = atom->x; double **f = atom->f; int nall = atom->nlocal + atom->nghost; @@ -237,8 +240,12 @@ void FixNumDiffVirial::calculate_virial() // restore original forces for owned and ghost atoms + atom->sync_host_arrays(F_MASK); + for (int i = 0; i < nall; i++) for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k]; + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- @@ -247,11 +254,19 @@ void FixNumDiffVirial::calculate_virial() void FixNumDiffVirial::displace_atoms(int nall, int idir, double magnitude) { + // the strain goes into the plain coordinate array, and the energy evaluation + // that follows works from the KOKKOS copies, so bring the host side up to + // date first and hand the write over afterwards + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int k = dirlist[idir][0]; int l = dirlist[idir][1]; for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k] + delta * magnitude * (temp_x[i][l] - fixedpoint[l]); + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -260,9 +275,13 @@ void FixNumDiffVirial::displace_atoms(int nall, int idir, double magnitude) void FixNumDiffVirial::restore_atoms(int nall, int idir) { + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int k = dirlist[idir][0]; for (int i = 0; i < nall; i++) { x[i][k] = temp_x[i][k]; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 886dfbec71a..c6204392750 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -577,6 +577,7 @@ action rand_pool_wrap_kokkos.cpp action rand_pool_wrap_kokkos.h action region_block_kokkos.cpp action region_block_kokkos.h +action region_remap_kokkos.h action region_sphere_kokkos.cpp action region_sphere_kokkos.h action remap_kokkos.cpp remap.cpp diff --git a/src/KOKKOS/atom_kokkos.cpp b/src/KOKKOS/atom_kokkos.cpp index ec737f4f6a7..cef9d75968f 100644 --- a/src/KOKKOS/atom_kokkos.cpp +++ b/src/KOKKOS/atom_kokkos.cpp @@ -275,6 +275,36 @@ void AtomKokkos::sync_pinned(const ExecutionSpace space, uint64_t mask, int asyn avecKK->sync_pinned(space, mask, async_flag); for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->sync_pinned(space, mask, async_flag); } +/* ---------------------------------------------------------------------- + the four ways of setting the per-type masses all write the plain host + array, so claim that write for the device copy +------------------------------------------------------------------------- */ + +void AtomKokkos::set_mass(const char *file, int line, const char *str, + int type_offset, int labelflag, int *ilabel) +{ + Atom::set_mass(file, line, str, type_offset, labelflag, ilabel); + k_mass.modify_host(); +} + +void AtomKokkos::set_mass(const char *file, int line, int itype, double value) +{ + Atom::set_mass(file, line, itype, value); + k_mass.modify_host(); +} + +void AtomKokkos::set_mass(const char *file, int line, int narg, char **arg) +{ + Atom::set_mass(file, line, narg, arg); + k_mass.modify_host(); +} + +void AtomKokkos::set_mass(double *values) +{ + Atom::set_mass(values); + k_mass.modify_host(); +} + /* ---------------------------------------------------------------------- */ void AtomKokkos::allocate_type_arrays() diff --git a/src/KOKKOS/atom_kokkos.h b/src/KOKKOS/atom_kokkos.h index 981b8dee91f..ddd0a508bec 100644 --- a/src/KOKKOS/atom_kokkos.h +++ b/src/KOKKOS/atom_kokkos.h @@ -189,6 +189,19 @@ class AtomKokkos : public Atom { void init() override; void update_property_atom(); void allocate_type_arrays() override; + + void sync_host_arrays(uint64_t mask) override { sync(Host, mask); } + void modified_host_arrays(uint64_t mask) override { modified(Host, mask); } + + // the per-type masses are written through the plain host array, which leaves + // the device copy behind with nothing to say so. Claim the write here, at + // the one place all four spellings of the mass command go through, rather + // than in each of the styles that read the masses on the device. + + void set_mass(const char *, int, const char *, int, int, int *) override; + void set_mass(const char *, int, int, double) override; + void set_mass(const char *, int, int, char **) override; + void set_mass(double *) override; void *extract(const char *) override; void sync(const ExecutionSpace space, uint64_t mask); void modified(const ExecutionSpace space, uint64_t mask); diff --git a/src/KOKKOS/atom_map_kokkos.cpp b/src/KOKKOS/atom_map_kokkos.cpp index 67cf8dfd218..fc636f4f6ca 100644 --- a/src/KOKKOS/atom_map_kokkos.cpp +++ b/src/KOKKOS/atom_map_kokkos.cpp @@ -426,6 +426,7 @@ void AtomKokkos::map_delete() { memoryKK->destroy_kokkos(k_sametag, sametag); sametag = nullptr; + max_same = 0; if (map_style == MAP_ARRAY) { memoryKK->destroy_kokkos(k_map_array, map_array); diff --git a/src/KOKKOS/atom_vec_angle_kokkos.cpp b/src/KOKKOS/atom_vec_angle_kokkos.cpp index 2b6ed3a3d1e..c0e9045d123 100644 --- a/src/KOKKOS/atom_vec_angle_kokkos.cpp +++ b/src/KOKKOS/atom_vec_angle_kokkos.cpp @@ -27,8 +27,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecAngleKokkos::AtomVecAngleKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecAngle(lmp) +AtomVecAngleKokkos::AtomVecAngleKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecAngle(lmp), molecule(nullptr), special(nullptr), + bond_atom(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), angle_atom3(nullptr) { } diff --git a/src/KOKKOS/atom_vec_bond_kokkos.cpp b/src/KOKKOS/atom_vec_bond_kokkos.cpp index 5e8a924c9e9..cdf6e46a1f9 100644 --- a/src/KOKKOS/atom_vec_bond_kokkos.cpp +++ b/src/KOKKOS/atom_vec_bond_kokkos.cpp @@ -26,8 +26,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecBondKokkos::AtomVecBondKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecBond(lmp) +AtomVecBondKokkos::AtomVecBondKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecBond(lmp), molecule(nullptr), special(nullptr), + bond_atom(nullptr) { } diff --git a/src/KOKKOS/atom_vec_dpd_kokkos.cpp b/src/KOKKOS/atom_vec_dpd_kokkos.cpp index be7a5d4f697..6b6995356b7 100644 --- a/src/KOKKOS/atom_vec_dpd_kokkos.cpp +++ b/src/KOKKOS/atom_vec_dpd_kokkos.cpp @@ -28,8 +28,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecDPDKokkos::AtomVecDPDKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecDPD(lmp) +AtomVecDPDKokkos::AtomVecDPDKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecDPD(lmp), duChem(nullptr) { } diff --git a/src/KOKKOS/atom_vec_full_kokkos.cpp b/src/KOKKOS/atom_vec_full_kokkos.cpp index e1e5de88f0d..6caaba7d007 100644 --- a/src/KOKKOS/atom_vec_full_kokkos.cpp +++ b/src/KOKKOS/atom_vec_full_kokkos.cpp @@ -26,8 +26,12 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecFullKokkos::AtomVecFullKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecFull(lmp) +AtomVecFullKokkos::AtomVecFullKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecFull(lmp), q(nullptr), molecule(nullptr), + special(nullptr), bond_atom(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), + angle_atom3(nullptr), dihedral_atom1(nullptr), dihedral_atom2(nullptr), + dihedral_atom3(nullptr), dihedral_atom4(nullptr), improper_atom1(nullptr), + improper_atom2(nullptr), improper_atom3(nullptr), improper_atom4(nullptr) { } diff --git a/src/KOKKOS/atom_vec_molecular_kokkos.cpp b/src/KOKKOS/atom_vec_molecular_kokkos.cpp index ae72ea7b164..5c74cb917f6 100644 --- a/src/KOKKOS/atom_vec_molecular_kokkos.cpp +++ b/src/KOKKOS/atom_vec_molecular_kokkos.cpp @@ -27,8 +27,12 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecMolecularKokkos::AtomVecMolecularKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecMolecular(lmp) +AtomVecMolecularKokkos::AtomVecMolecularKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecMolecular(lmp), molecule(nullptr), special(nullptr), + bond_atom(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), angle_atom3(nullptr), + dihedral_atom1(nullptr), dihedral_atom2(nullptr), dihedral_atom3(nullptr), + dihedral_atom4(nullptr), improper_atom1(nullptr), improper_atom2(nullptr), + improper_atom3(nullptr), improper_atom4(nullptr) { } diff --git a/src/KOKKOS/atom_vec_sphere_kokkos.cpp b/src/KOKKOS/atom_vec_sphere_kokkos.cpp index 83a8ff37ff4..a40a246daec 100644 --- a/src/KOKKOS/atom_vec_sphere_kokkos.cpp +++ b/src/KOKKOS/atom_vec_sphere_kokkos.cpp @@ -30,8 +30,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -AtomVecSphereKokkos::AtomVecSphereKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecSphere(lmp) +AtomVecSphereKokkos::AtomVecSphereKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecSphere(lmp), torque(nullptr) { } diff --git a/src/KOKKOS/bond_quartic_kokkos.cpp b/src/KOKKOS/bond_quartic_kokkos.cpp index 38d12ae402c..14058d8a142 100644 --- a/src/KOKKOS/bond_quartic_kokkos.cpp +++ b/src/KOKKOS/bond_quartic_kokkos.cpp @@ -19,6 +19,7 @@ #include "bond_quartic_kokkos.h" #include "atom_kokkos.h" +#include "kokkos.h" #include "atom_masks.h" #include "comm.h" #include "force.h" @@ -44,6 +45,15 @@ BondQuarticKokkos::BondQuarticKokkos(LAMMPS *lmp) : BondQuartic(lmp) execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; + + // this style breaks bonds and corrects the 1-4 pair interaction on the host, + // in the middle of its own compute(), which means copying the forces to the + // host and back again. The overlap path keeps a separate host force buffer, + // zeroed and filled by the host executing styles alone and merged into the + // device buffer afterwards, and a copy in either direction from here lands + // in the middle of that. Turn the overlap off while this style is in use, + // as bond/hybrid does for a host executing sub-style. + lmp->kokkos->allow_overlap = 0; } /* ---------------------------------------------------------------------- */ @@ -90,6 +100,13 @@ void BondQuarticKokkos::compute(int eflag_in, int vflag_in) k_rc.template sync(); k_u0.template sync(); + // the kernel adds to the forces and the host pass below then copies them + // away and back, so they have to be current on this side before it runs + // whatever the caller did; on any path that reaches compute() without the + // integrator's own sync, claiming the device side after the kernel would + // otherwise collide with a pending host modification + atomKK->sync(execution_space,X_MASK|F_MASK); + x = atomKK->k_x.template view(); f = atomKK->k_f.template view(); neighborKK->k_bondlist.template sync(); diff --git a/src/KOKKOS/comm_tiled_kokkos.cpp b/src/KOKKOS/comm_tiled_kokkos.cpp index abfc1884bfe..62ecde2b7fa 100644 --- a/src/KOKKOS/comm_tiled_kokkos.cpp +++ b/src/KOKKOS/comm_tiled_kokkos.cpp @@ -668,6 +668,12 @@ void CommTiledKokkos::grow_buf_pair(int n) max_buf_pair = n * BUFFACTOR; k_buf_send_pair.resize(max_buf_pair); k_buf_recv_pair.resize(max_buf_pair); + + // resizing claims a side; these are scratch buffers that are filled + // before they are read, so drop the claim rather than leave it for the + // next modify_host() to collide with + k_buf_send_pair.clear_sync_state(); + k_buf_recv_pair.clear_sync_state(); } /* ---------------------------------------------------------------------- @@ -855,6 +861,12 @@ void CommTiledKokkos::grow_send_kokkos(int n, int flag, ExecutionSpace space) atomKK->avecKK->size_border + atomKK->avecKK->size_velocity); else k_buf_send.resize(maxsend_border,atomKK->avecKK->size_border); + + // the claim above only steers the resize to the side whose contents have + // to survive; after it this is a scratch buffer again, filled through raw + // pointers on whichever side does the packing, so drop the claim rather + // than leave it standing forever + k_buf_send.clear_sync_state(); } else { if (ghost_velocity) MemoryKokkos::realloc_kokkos(k_buf_send,"comm:k_buf_send",maxsend_border, diff --git a/src/KOKKOS/domain_kokkos.cpp b/src/KOKKOS/domain_kokkos.cpp index 1d8718c0694..a2e1caaba6a 100644 --- a/src/KOKKOS/domain_kokkos.cpp +++ b/src/KOKKOS/domain_kokkos.cpp @@ -416,6 +416,17 @@ void DomainKokkos::pbc() void DomainKokkos::remap_all() { + // the kernel below does not carry the velocity correction that fix deform's + // "remap v" applies when an atom is wrapped, so leave that case to the host + // implementation rather than drop the correction + + if (deform_vremap) { + atomKK->sync(Host,X_MASK | V_MASK | IMAGE_MASK | MASK_MASK); + Domain::remap_all(); + atomKK->modified(Host,X_MASK | V_MASK | IMAGE_MASK); + return; + } + atomKK->sync(Device,X_MASK | IMAGE_MASK); x = atomKK->k_x.view_device(); diff --git a/src/KOKKOS/domain_kokkos.h b/src/KOKKOS/domain_kokkos.h index e94e3f13122..2c5775f96e8 100644 --- a/src/KOKKOS/domain_kokkos.h +++ b/src/KOKKOS/domain_kokkos.h @@ -34,8 +34,8 @@ class DomainKokkos : public Domain { ~DomainKokkos() override = default; void reset_box() override; void pbc() override; - void remap_all(); - void image_flip(int, int, int); + void remap_all() override; + void image_flip(int, int, int) override; void x2lamda(int) override; void x2lamda(int,int) override; void lamda2x(int) override; diff --git a/src/KOKKOS/fix_addforce_kokkos.cpp b/src/KOKKOS/fix_addforce_kokkos.cpp index 101c039be6d..446bc015ed0 100644 --- a/src/KOKKOS/fix_addforce_kokkos.cpp +++ b/src/KOKKOS/fix_addforce_kokkos.cpp @@ -74,6 +74,8 @@ void FixAddForceKokkos::init() template void FixAddForceKokkos::post_force(int vflag) { + if (update->ntimestep % nevery) return; + atomKK->sync(execution_space, X_MASK | F_MASK | IMAGE_MASK | MASK_MASK); x = atomKK->k_x.view(); @@ -84,13 +86,16 @@ void FixAddForceKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); - memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"efield:vatom"); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"addforce:vatom"); d_vatom = k_vatom.template view(); } @@ -109,7 +114,7 @@ void FixAddForceKokkos::post_force(int vflag) // reallocate sforce array if necessary - if (varflag == ATOM && atom->nmax > maxatom) { + if (((varflag == ATOM) || (estyle == ATOM)) && atom->nmax > maxatom) { maxatom = atom->nmax; memoryKK->destroy_kokkos(k_sforce,sforce); memoryKK->create_kokkos(k_sforce,sforce,maxatom,4,"addforce:sforce"); @@ -150,8 +155,10 @@ void FixAddForceKokkos::post_force(int vflag) // atom-style variables are evaluated on the host, so the result has to be // copied to the device for the kernel below. this is a real copy, not a // stale flag: it can only go away if variables are evaluated on the device. + // an atom-style energy variable writes into sforce[i][3], so it needs the + // copy as well - if (varflag == ATOM) { + if ((varflag == ATOM) || (estyle == ATOM)) { k_sforce.modify_host(); k_sforce.sync(); } @@ -206,6 +213,17 @@ void FixAddForceKokkos::operator()(TagFixAddForceConstant, const int if (xstyle) f(i,0) += static_cast(xvalue_kk); if (ystyle) f(i,1) += static_cast(yvalue_kk); if (zstyle) f(i,2) += static_cast(zvalue_kk); + + if (evflag) { + KK_FLOAT v[6]; + v[0] = static_cast(xvalue * unwrapKK[0]); + v[1] = static_cast(yvalue * unwrapKK[1]); + v[2] = static_cast(zvalue * unwrapKK[2]); + v[3] = static_cast(xvalue * unwrapKK[1]); + v[4] = static_cast(xvalue * unwrapKK[2]); + v[5] = static_cast(yvalue * unwrapKK[2]); + v_tally(result,i,v); + } } } @@ -225,15 +243,19 @@ void FixAddForceKokkos::operator()(TagFixAddForceNonConstant, const x_i[2] = static_cast(x(i,2)); auto unwrapKK = DomainKokkos::unmap(prd,h,triclinic,x_i,image(i)); + // an atom-style variable supplies the value per atom, any other style + // (constant or equal-style variable) the same value for all of them + + const double xv = (xstyle == ATOM) ? static_cast(d_sforce(i,0)) : xvalue; + const double yv = (ystyle == ATOM) ? static_cast(d_sforce(i,1)) : yvalue; + const double zv = (zstyle == ATOM) ? static_cast(d_sforce(i,2)) : zvalue; + if (estyle == ATOM) { result[0] += static_cast(d_sforce(i,3)); } else { - if (xstyle == EQUAL) result[0] -= xvalue * unwrapKK[0]; - if (ystyle == EQUAL) result[0] -= yvalue * unwrapKK[1]; - if (zstyle == EQUAL) result[0] -= zvalue * unwrapKK[2]; - if (xstyle == ATOM) result[0] -= static_cast(d_sforce(i,0)) * unwrapKK[0]; - if (ystyle == ATOM) result[0] -= static_cast(d_sforce(i,1)) * unwrapKK[1]; - if (zstyle == ATOM) result[0] -= static_cast(d_sforce(i,2)) * unwrapKK[2]; + if (xstyle) result[0] -= xv * unwrapKK[0]; + if (ystyle) result[0] -= yv * unwrapKK[1]; + if (zstyle) result[0] -= zv * unwrapKK[2]; } result[1] += static_cast(f(i,0)); result[2] += static_cast(f(i,1)); @@ -244,6 +266,17 @@ void FixAddForceKokkos::operator()(TagFixAddForceNonConstant, const else if (ystyle) f(i,1) += static_cast(yvalue_kk); if (zstyle == ATOM) f(i,2) += static_cast(d_sforce(i,2)); else if (zstyle) f(i,2) += static_cast(zvalue_kk); + + if (evflag) { + KK_FLOAT v[6]; + v[0] = xstyle ? static_cast(xv * unwrapKK[0]) : static_cast(0.0); + v[1] = ystyle ? static_cast(yv * unwrapKK[1]) : static_cast(0.0); + v[2] = zstyle ? static_cast(zv * unwrapKK[2]) : static_cast(0.0); + v[3] = xstyle ? static_cast(xv * unwrapKK[1]) : static_cast(0.0); + v[4] = xstyle ? static_cast(xv * unwrapKK[2]) : static_cast(0.0); + v[5] = ystyle ? static_cast(yv * unwrapKK[2]) : static_cast(0.0); + v_tally(result,i,v); + } } } diff --git a/src/KOKKOS/fix_deform_kokkos.cpp b/src/KOKKOS/fix_deform_kokkos.cpp index 5494748f3a3..b0c5c271dd8 100644 --- a/src/KOKKOS/fix_deform_kokkos.cpp +++ b/src/KOKKOS/fix_deform_kokkos.cpp @@ -46,8 +46,23 @@ FixDeformKokkos::FixDeformKokkos(LAMMPS *lmp, int narg, char **arg) : FixDeform( void FixDeformKokkos::pre_exchange() { - atomKK->sync(Host,ALL_MASK); + // everything the base class does here goes through DomainKokkos, which runs + // on the device and declares itself, except the atom migration below. + // Bracketing the whole call for the host claimed a side the device work had + // already claimed, and left the migration reading host data that the device + // work had since overtaken. + FixDeform::pre_exchange(); +} + +/* ---------------------------------------------------------------------- + the migration is the one host only step in pre_exchange() +------------------------------------------------------------------------- */ + +void FixDeformKokkos::migrate_atoms() +{ + atomKK->sync(Host,ALL_MASK); + FixDeform::migrate_atoms(); atomKK->modified(Host,ALL_MASK); } diff --git a/src/KOKKOS/fix_deform_kokkos.h b/src/KOKKOS/fix_deform_kokkos.h index a306a291ef4..0c30f85ffe2 100644 --- a/src/KOKKOS/fix_deform_kokkos.h +++ b/src/KOKKOS/fix_deform_kokkos.h @@ -32,6 +32,7 @@ class FixDeformKokkos : public FixDeform { FixDeformKokkos(class LAMMPS *, int, char **); void pre_exchange() override; + void migrate_atoms() override; void update_box() override; }; diff --git a/src/KOKKOS/fix_efield_kokkos.cpp b/src/KOKKOS/fix_efield_kokkos.cpp index 00a1a14f162..f3746cd68f4 100644 --- a/src/KOKKOS/fix_efield_kokkos.cpp +++ b/src/KOKKOS/fix_efield_kokkos.cpp @@ -93,9 +93,12 @@ void FixEfieldKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) - // reallocate per-atom arrays if necessary + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); @@ -116,9 +119,11 @@ void FixEfieldKokkos::post_force(int vflag) d_match = k_match.template view(); } - // reallocate sforce array if necessary + // reallocate efield array if necessary + // an atom-style energy or potential variable writes into efield[i][3], + // so the array is needed for those as well, not only for varflag == ATOM - if (varflag == ATOM && atom->nmax > maxatom) { + if (((varflag == ATOM) || (estyle == ATOM) || (pstyle == ATOM)) && atom->nmax > maxatom) { maxatom = atom->nmax; memoryKK->destroy_kokkos(k_efield,efield); memoryKK->create_kokkos(k_efield,efield,maxatom,4,"efield:efield"); @@ -155,8 +160,10 @@ void FixEfieldKokkos::post_force(int vflag) // atom-style variables are evaluated on the host, so the result has to be // copied to the device for the kernel below. this is a real copy, not a // stale flag: it can only go away if variables are evaluated on the device. + // an atom-style energy or potential variable writes into efield[i][3], so + // it needs the copy as well - if (varflag == ATOM) { + if ((varflag == ATOM) || (estyle == ATOM) || (pstyle == ATOM)) { k_efield.modify_host(); k_efield.sync(); } diff --git a/src/KOKKOS/fix_minimize_kokkos.cpp b/src/KOKKOS/fix_minimize_kokkos.cpp index b1d423c89d7..741e46e7fc9 100644 --- a/src/KOKKOS/fix_minimize_kokkos.cpp +++ b/src/KOKKOS/fix_minimize_kokkos.cpp @@ -31,6 +31,12 @@ FixMinimizeKokkos::FixMinimizeKokkos(LAMMPS *lmp, int narg, char **arg) : { kokkosable = 1; atomKK = (AtomKokkos *) atom; + + // this fix only stores per-atom data of its own; it syncs the atom data it + // touches itself, so Modify must not sync or invalidate any of it + + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/fix_nh_kokkos.cpp b/src/KOKKOS/fix_nh_kokkos.cpp index 0e8948de7ce..0ad9712b2f9 100644 --- a/src/KOKKOS/fix_nh_kokkos.cpp +++ b/src/KOKKOS/fix_nh_kokkos.cpp @@ -325,13 +325,17 @@ void FixNHKokkos::remap() if (allremap) domainKK->x2lamda(nlocal); else { - for ( int i = 0; i < nlocal; i++) - if (mask[i] & dilate_group_bit) { - auto h_x = atomKK->k_x.view_host(); - atomKK->sync(Host,X_MASK); + // this loop runs on the host, so it needs the host side of both arrays: + // "mask" is the view for the execution space, and reading it here is an + // access to device memory from host code. Sync once for the whole loop + // rather than once per atom as well. + atomKK->sync(Host,X_MASK|MASK_MASK); + auto h_x = atomKK->k_x.view_host(); + auto h_mask = atomKK->k_mask.view_host(); + for (int i = 0; i < nlocal; i++) + if (h_mask[i] & dilate_group_bit) domainKK->x2lamda(&h_x(i,0), &h_x(i,0)); - atomKK->modified(Host,X_MASK); - } + atomKK->modified(Host,X_MASK); } if (rfix.size() > 0) @@ -476,13 +480,17 @@ void FixNHKokkos::remap() if (allremap) domainKK->lamda2x(nlocal); else { - for ( int i = 0; i < nlocal; i++) - if (mask[i] & dilate_group_bit) { - auto h_x = atomKK->k_x.view_host(); - atomKK->sync(Host,X_MASK); + // this loop runs on the host, so it needs the host side of both arrays: + // "mask" is the view for the execution space, and reading it here is an + // access to device memory from host code. Sync once for the whole loop + // rather than once per atom as well. + atomKK->sync(Host,X_MASK|MASK_MASK); + auto h_x = atomKK->k_x.view_host(); + auto h_mask = atomKK->k_mask.view_host(); + for (int i = 0; i < nlocal; i++) + if (h_mask[i] & dilate_group_bit) domainKK->lamda2x(&h_x(i,0), &h_x(i,0)); - atomKK->modified(Host,X_MASK); - } + atomKK->modified(Host,X_MASK); } // for (auto &ifix : rfix) ifix->deform(1); @@ -618,7 +626,6 @@ template void FixNHKokkos::nve_x() { atomKK->sync(execution_space,X_MASK | V_MASK | MASK_MASK); - atomKK->modified(execution_space,X_MASK); x = atomKK->k_x.view(); v = atomKK->k_v.view(); @@ -631,6 +638,12 @@ void FixNHKokkos::nve_x() copymode = 1; Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); copymode = 0; + + // claim the coordinates after the kernel has written them, not before: a + // claim made up front can be taken by a copy that runs in between, which + // leaves the coordinates this kernel writes unclaimed + + atomKK->modified(execution_space,X_MASK); } template diff --git a/src/KOKKOS/fix_property_atom_kokkos.cpp b/src/KOKKOS/fix_property_atom_kokkos.cpp index ecb60495b14..fc003b71711 100644 --- a/src/KOKKOS/fix_property_atom_kokkos.cpp +++ b/src/KOKKOS/fix_property_atom_kokkos.cpp @@ -31,6 +31,12 @@ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) atomKK = (AtomKokkos *) atom; kokkosable = 1; + // this fix syncs the atom data it owns itself, through AtomKokkos::sync() + // and AtomKokkos::modified(); Modify must not sync or invalidate any of it + + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; + dvector_flag = 0; ivector_flag = 0; iarray_flag = 0; diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 4fff360a394..ca5c10eb523 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -49,6 +49,12 @@ FixShakeKokkos::FixShakeKokkos(LAMMPS *lmp, int narg, char **arg) : if (store_flag) error->all(FLERR, "Option 'store yes' is not (yet) supported by fix {}/kk", style); + // the centroid virial of the SHAKE constraint forces is not (yet) computed + // by the Kokkos kernels, so do not claim it is available. compute + // centroid/stress/atom then errors out instead of silently reporting zeros. + + centroidstressflag = CENTROID_NOTAVAIL; + datamask_read = EMPTY_MASK; datamask_modify = EMPTY_MASK; @@ -141,6 +147,7 @@ FixShakeKokkos::~FixShakeKokkos() memoryKK->destroy_kokkos(k_list,list); memoryKK->destroy_kokkos(k_closest_list,closest_list); + memoryKK->destroy_kokkos(k_eatom,eatom); memoryKK->destroy_kokkos(k_vatom,vatom); } @@ -152,6 +159,15 @@ FixShakeKokkos::~FixShakeKokkos() template void FixShakeKokkos::init() { + // FixShake::init() reads shake_flag/shake_atom/shake_type through the + // plain host pointers to recompute angle_distance[]; comm->exchange() + // updates the device copies only, so the host side must be brought + // current before FixShake::init() can see this run's cluster data + + k_shake_flag.sync_host(); + k_shake_atom.sync_host(); + k_shake_type.sync_host(); + FixShake::init(); if (utils::strmatch(update->integrate_style,"^respa")) @@ -345,17 +361,49 @@ template void FixShakeKokkos::min_post_force(int vflag) { int eflag = eflag_pre_reverse; - ev_init(eflag, vflag); + + // the per-atom virial is accumulated into a dual view, so the plain + // base-class arrays must not be allocated here (alloc = 0) + + ev_init(eflag, vflag, 0); + + // reallocate the per-atom energy and virial dual views if necessary + + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"shake:eatom"); + d_eatom = k_eatom.template view(); + } + if (vflag_atom) { + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"shake:vatom"); + d_vatom = k_vatom.template view(); + } + ebond = 0.0; + nlocal = atomKK->nlocal; atomKK->sync(execution_space, X_MASK | F_MASK); k_shake_flag.sync(); + k_shake_atom.sync(); k_shake_type.sync(); k_list.sync(); k_closest_list.sync(); k_bond_distance.sync(); k_angle_distance.sync(); + // the restraint energy and virial are credited to the owner of each atom, + // which requires mapping image atoms back through the atom map + + map_style = atom->map_style; + if (map_style == Atom::MAP_ARRAY) { + k_map_array = atomKK->k_map_array; + k_map_array.template sync(); + } else if (map_style == Atom::MAP_HASH) { + k_map_hash = atomKK->k_map_hash; + k_map_hash.template sync(); + } + // Assign class member views for Tagged Kernel access this->d_x = atomKK->k_x.view(); this->d_f = atomKK->k_f.view(); @@ -369,10 +417,15 @@ void FixShakeKokkos::min_post_force(int vflag) if (neighflag != HALF) need_dup = std::is_same_v, Kokkos::Experimental::ScatterDuplicated>; - if (need_dup) - dup_f = Kokkos::Experimental::create_scatter_view(d_f); - else - ndup_f = Kokkos::Experimental::create_scatter_view(d_f); + if (need_dup) { + dup_f = Kokkos::Experimental::create_scatter_view(d_f); + dup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + dup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } else { + ndup_f = Kokkos::Experimental::create_scatter_view(d_f); + ndup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + ndup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } copymode = 1; @@ -396,20 +449,27 @@ void FixShakeKokkos::min_post_force(int vflag) EV_FLOAT ev; if (neighflag == HALF) { - if(vflag) + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); else - Kokkos::parallel_for(Kokkos::RangePolicy>(0, nlist), *this); + Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); } else { - if(vflag) + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); else - Kokkos::parallel_for(Kokkos::RangePolicy>(0, nlist), *this); + Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); } copymode = 0; - if (need_dup) Kokkos::Experimental::contribute(d_f, dup_f); + // reduction over duplicated memory + + if (need_dup) { + Kokkos::Experimental::contribute(d_f, dup_f); + if (eflag_atom) Kokkos::Experimental::contribute(d_eatom, dup_eatom); + if (vflag_atom) Kokkos::Experimental::contribute(d_vatom, dup_vatom); + } + comm->reverse_comm(this); this->ebond = static_cast(ev.evdwl); @@ -422,22 +482,37 @@ void FixShakeKokkos::min_post_force(int vflag) virial[5] += static_cast(ev.v[5]); } + if (eflag_atom) { + k_eatom.template modify(); + k_eatom.sync_host(); + } + if (vflag_atom) { + k_vatom.template modify(); + k_vatom.sync_host(); + } + atomKK->modified(execution_space, F_MASK); + // free duplicated memory + + if (need_dup) { + dup_f = {}; + dup_eatom = {}; + dup_vatom = {}; + } + if (update->ntimestep == next_output) { atomKK->modified(execution_space, X_MASK); stats(); } - - if (need_dup) dup_f = {}; } /* ---------------------------------------------------------------------- */ template -template +template KOKKOS_INLINE_FUNCTION -void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i, EV_FLOAT &ev) const +void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i, EV_FLOAT &ev) const { auto v_f = ScatterViewHelper, decltype(dup_f), decltype(ndup_f)>::get(dup_f, ndup_f); auto a_f = v_f.template access>(); @@ -446,7 +521,12 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce::operator()(TagFixShakeMinPostForce(delx * fbond); a_f(idx1, 1) -= static_cast(dely * fbond); a_f(idx1, 2) -= static_cast(delz * fbond); - ev.evdwl += static_cast(eb); - if (VFLAG) { - ev.v[0] += static_cast(static_cast(0.5) * delx * delx * fbond); - ev.v[1] += static_cast(static_cast(0.5) * dely * dely * fbond); - ev.v[2] += static_cast(static_cast(0.5) * delz * delz * fbond); - ev.v[3] += static_cast(static_cast(0.5) * delx * dely * fbond); - ev.v[4] += static_cast(static_cast(0.5) * delx * delz * fbond); - ev.v[5] += static_cast(static_cast(0.5) * dely * delz * fbond); + + // energy and virial are shared out over the owned atoms of the pair, as + // FixShake::bond_force() does. the closest image of an atom owned by + // this processor is mapped back onto the owner first, so that a cluster + // reaching across a periodic boundary is still credited to both atoms. + + int atomlist[2]; + int count = 0; + const int own0 = AtomKokkos::map_kokkos(d_shake_atom(m,slot0),map_style, + k_map_array,k_map_hash); + const int own1 = AtomKokkos::map_kokkos(d_shake_atom(m,slot1),map_style, + k_map_array,k_map_hash); + if ((own0 >= 0) && (own0 < nlocal)) atomlist[count++] = own0; + if ((own1 >= 0) && (own1 < nlocal)) atomlist[count++] = own1; + + const KK_FLOAT total = static_cast(2.0); + ev.evdwl += static_cast((static_cast(count)/total) * eb); + + if (EVFLAG) { + KK_FLOAT v[6]; + v[0] = static_cast(0.5) * delx * delx * fbond; + v[1] = static_cast(0.5) * dely * dely * fbond; + v[2] = static_cast(0.5) * delz * delz * fbond; + v[3] = static_cast(0.5) * delx * dely * fbond; + v[4] = static_cast(0.5) * delx * delz * fbond; + v[5] = static_cast(0.5) * dely * delz * fbond; + ev_tally(ev,count,atomlist,total,eb,v); } if (output_every && !is_angle) { - Kokkos::atomic_add(&d_b_stats(type_idx, 0), 1.0); - Kokkos::atomic_add(&d_b_stats(type_idx, 1), (double)r); - Kokkos::atomic_add(&d_b_stats(type_idx, 0), 1.0); - Kokkos::atomic_add(&d_b_stats(type_idx, 1), (double)r); + // only atoms owned by this processor are counted, as on the CPU + const double nown = (double) ((idx0 < nlocal) + (idx1 < nlocal)); + if (nown > 0.0) { + Kokkos::atomic_add(&d_b_stats(type_idx, 0), nown); + Kokkos::atomic_add(&d_b_stats(type_idx, 1), nown * (double)r); + } Kokkos::atomic_max(&d_b_stats(type_idx, 2), (double)r); Kokkos::atomic_min(&d_b_stats(type_idx, 3), (double)r); } @@ -484,20 +585,20 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce(2.0)*r1*r2)) * static_cast(180.0)/static_cast(MY_PI); int mt = d_shake_type(m, 2); @@ -513,11 +614,11 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce -template +template KOKKOS_INLINE_FUNCTION -void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i) const { +void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i) const { EV_FLOAT ev; - this->template operator()(TagFixShakeMinPostForce(), i, ev); + this->template operator()(TagFixShakeMinPostForce(), i, ev); } /* ---------------------------------------------------------------------- @@ -566,15 +667,26 @@ void FixShakeKokkos::post_force(int vflag) comm->forward_comm(this); k_xshake.sync(); - // virial setup + // energy and virial setup, as FixShake::post_force() does. the per-atom + // energy and virial are accumulated into dual views, so the plain + // base-class arrays must not be allocated here (alloc = 0) - v_init(vflag); + int eflag = eflag_pre_reverse; + ev_init(eflag,vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom energy and virial dual views if necessary. the + // constraint forces contribute no per-atom energy during dynamics, but the + // freshly created view keeps it zeroed rather than holding on to whatever + // a preceding minimization left behind + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"shake:eatom"); + d_eatom = k_eatom.template view(); + } if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); - memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"improper:vatom"); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"shake:vatom"); d_vatom = k_vatom.template view(); } @@ -639,8 +751,10 @@ void FixShakeKokkos::post_force(int vflag) // reduction over duplicated memory - if (need_dup) + if (need_dup) { Kokkos::Experimental::contribute(d_f,dup_f); + if (vflag_atom) Kokkos::Experimental::contribute(d_vatom,dup_vatom); + } atomKK->modified(execution_space,F_MASK); @@ -1617,9 +1731,13 @@ void FixShakeKokkos::stats() template void FixShakeKokkos::grow_arrays(int nmax) { - k_shake_flag.sync_device(); - k_shake_atom.sync_device(); - k_shake_type.sync_device(); + // grow_kokkos() grows the host side and hands the plain pointers below it + // back to the base class, which reads and writes the clusters through them, + // so the host side is the one that has to be current here + + k_shake_flag.sync_host(); + k_shake_atom.sync_host(); + k_shake_type.sync_host(); memoryKK->grow_kokkos(k_shake_flag,shake_flag,nmax,"shake:shake_flag"); memoryKK->grow_kokkos(k_shake_atom,shake_atom,nmax,4,"shake:shake_atom"); @@ -1988,6 +2106,13 @@ int FixShakeKokkos::pack_forward_comm_kokkos(int n, DAT::tdual_int_1 DAT::tdual_double_1d &k_buf, int pbc_flag, int* pbc) { + // the host variant below syncs before it reads and claims after it writes; + // this one has to do the same, or the buffer is packed from whichever side + // the device copy last held and the ghosts it unpacks are discarded by the + // next sync + + k_xshake.sync(); + d_sendlist = k_sendlist.view(); d_buf = k_buf.view(); @@ -2049,6 +2174,8 @@ void FixShakeKokkos::unpack_forward_comm_kokkos(int n, int first_in, first = first_in; d_buf = buf.view(); Kokkos::parallel_for(Kokkos::RangePolicy(0,n),*this); + + k_xshake.modify(); } template @@ -2183,6 +2310,32 @@ void FixShakeKokkos::correct_coordinates(int vflag) { atomKK->modified(Host,X_MASK|V_MASK|F_MASK); } +/* ---------------------------------------------------------------------- + tally energy and virial into global and per-atom accumulators + n = # of local owned atoms involved, with local indices in list + eng = total energy for the interaction involving total atoms + increment per-atom energy of each atom in list by 1/total fraction + mirrors Fix::ev_tally() +------------------------------------------------------------------------- */ + +template +template +// NOLINTNEXTLINE +KOKKOS_INLINE_FUNCTION +void FixShakeKokkos::ev_tally(EV_FLOAT &ev, int n, int *atomlist, KK_FLOAT total, + KK_FLOAT eng, KK_FLOAT *v) const +{ + if (eflag_atom) { + auto v_eatom = ScatterViewHelper,decltype(dup_eatom),decltype(ndup_eatom)>::get(dup_eatom,ndup_eatom); + auto a_eatom = v_eatom.template access>(); + const KK_FLOAT fraction = eng/total; + for (int i = 0; i < n; i++) + a_eatom(atomlist[i]) += static_cast(fraction); + } + + v_tally(ev,n,atomlist,total,v); +} + /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators n = # of local owned atoms involved, with local indices in list diff --git a/src/KOKKOS/fix_shake_kokkos.h b/src/KOKKOS/fix_shake_kokkos.h index 13a4dd6f896..afea28d331f 100644 --- a/src/KOKKOS/fix_shake_kokkos.h +++ b/src/KOKKOS/fix_shake_kokkos.h @@ -32,7 +32,7 @@ namespace LAMMPS_NS { struct TagFixShakePreNeighbor{}; -template +template struct TagFixShakeMinPostForce{}; template @@ -86,13 +86,13 @@ class FixShakeKokkos : public FixShake, public KokkosBase { KOKKOS_INLINE_FUNCTION void operator()(TagFixShakePreNeighbor, const int&) const; - template + template KOKKOS_INLINE_FUNCTION - void operator()(TagFixShakeMinPostForce, const int&, EV_FLOAT&) const; + void operator()(TagFixShakeMinPostForce, const int&, EV_FLOAT&) const; - template + template KOKKOS_INLINE_FUNCTION - void operator()(TagFixShakeMinPostForce, const int&) const; + void operator()(TagFixShakeMinPostForce, const int&) const; template // NOLINTNEXTLINE @@ -234,6 +234,11 @@ class FixShakeKokkos : public FixShake, public KokkosBase { KOKKOS_INLINE_FUNCTION void v_tally(EV_FLOAT&, int, int *, KK_FLOAT, KK_FLOAT *) const; + template +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void ev_tally(EV_FLOAT&, int, int *, KK_FLOAT, KK_FLOAT, KK_FLOAT *) const; + int first,nsend; typename AT::t_int_1d d_sendlist; diff --git a/src/KOKKOS/fix_wall_harmonic_kokkos.cpp b/src/KOKKOS/fix_wall_harmonic_kokkos.cpp index 49b9dc299f5..fb3c1a90b29 100644 --- a/src/KOKKOS/fix_wall_harmonic_kokkos.cpp +++ b/src/KOKKOS/fix_wall_harmonic_kokkos.cpp @@ -55,14 +55,27 @@ void FixWallHarmonicKokkos::precompute(int /*m_in*/) /* ---------------------------------------------------------------------- */ template -void FixWallHarmonicKokkos::post_force(int vflag) +void FixWallHarmonicKokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_harmonic:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallHarmonicKokkos::post_force(int vflag) +{ FixWallHarmonic::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_harmonic_kokkos.h b/src/KOKKOS/fix_wall_harmonic_kokkos.h index 582a4b664a4..97d0f4955ff 100644 --- a/src/KOKKOS/fix_wall_harmonic_kokkos.h +++ b/src/KOKKOS/fix_wall_harmonic_kokkos.h @@ -40,6 +40,7 @@ class FixWallHarmonicKokkos : public FixWallHarmonic { ~FixWallHarmonicKokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_lj1043_kokkos.cpp b/src/KOKKOS/fix_wall_lj1043_kokkos.cpp index ba2b5c34e4e..d18baf32028 100644 --- a/src/KOKKOS/fix_wall_lj1043_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj1043_kokkos.cpp @@ -120,14 +120,27 @@ void FixWallLJ1043Kokkos::precompute(int m_in) /* ---------------------------------------------------------------------- */ template -void FixWallLJ1043Kokkos::post_force(int vflag) +void FixWallLJ1043Kokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_lj1043:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallLJ1043Kokkos::post_force(int vflag) +{ FixWallLJ1043::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_lj1043_kokkos.h b/src/KOKKOS/fix_wall_lj1043_kokkos.h index 31ed5e97a49..d254e86b31e 100644 --- a/src/KOKKOS/fix_wall_lj1043_kokkos.h +++ b/src/KOKKOS/fix_wall_lj1043_kokkos.h @@ -40,6 +40,7 @@ class FixWallLJ1043Kokkos : public FixWallLJ1043 { ~FixWallLJ1043Kokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_lj126_kokkos.cpp b/src/KOKKOS/fix_wall_lj126_kokkos.cpp index 1ced7256191..a2704bc4010 100644 --- a/src/KOKKOS/fix_wall_lj126_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj126_kokkos.cpp @@ -98,14 +98,27 @@ void FixWallLJ126Kokkos::precompute(int m_in) /* ---------------------------------------------------------------------- */ template -void FixWallLJ126Kokkos::post_force(int vflag) +void FixWallLJ126Kokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_lj126:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallLJ126Kokkos::post_force(int vflag) +{ FixWallLJ126::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_lj126_kokkos.h b/src/KOKKOS/fix_wall_lj126_kokkos.h index e86dca95065..6120790419d 100644 --- a/src/KOKKOS/fix_wall_lj126_kokkos.h +++ b/src/KOKKOS/fix_wall_lj126_kokkos.h @@ -40,6 +40,7 @@ class FixWallLJ126Kokkos : public FixWallLJ126 { ~FixWallLJ126Kokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.cpp b/src/KOKKOS/fix_wall_lj93_kokkos.cpp index d63ebfc4caf..f7c4c6f8539 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj93_kokkos.cpp @@ -105,17 +105,27 @@ void FixWallLJ93Kokkos::precompute(int m) /* ---------------------------------------------------------------------- */ template -void FixWallLJ93Kokkos::post_force(int vflag) +void FixWallLJ93Kokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"wall_lj93:vatom"); d_vatom = k_vatom.template view(); } +} + +/* ---------------------------------------------------------------------- */ +template +void FixWallLJ93Kokkos::post_force(int vflag) +{ FixWallLJ93::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.h b/src/KOKKOS/fix_wall_lj93_kokkos.h index 51bdd7dbd2c..334348490a5 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.h +++ b/src/KOKKOS/fix_wall_lj93_kokkos.h @@ -40,6 +40,7 @@ class FixWallLJ93Kokkos : public FixWallLJ93 { ~FixWallLJ93Kokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_morse_kokkos.cpp b/src/KOKKOS/fix_wall_morse_kokkos.cpp index 1f547f381f0..396422eff4a 100644 --- a/src/KOKKOS/fix_wall_morse_kokkos.cpp +++ b/src/KOKKOS/fix_wall_morse_kokkos.cpp @@ -100,14 +100,27 @@ void FixWallMorseKokkos::precompute(int m_in) /* ---------------------------------------------------------------------- */ template -void FixWallMorseKokkos::post_force(int vflag) +void FixWallMorseKokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_morse:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallMorseKokkos::post_force(int vflag) +{ FixWallMorse::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_morse_kokkos.h b/src/KOKKOS/fix_wall_morse_kokkos.h index 3495549e35d..f1cf9197578 100644 --- a/src/KOKKOS/fix_wall_morse_kokkos.h +++ b/src/KOKKOS/fix_wall_morse_kokkos.h @@ -40,6 +40,7 @@ class FixWallMorseKokkos : public FixWallMorse { ~FixWallMorseKokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_region_kokkos.cpp b/src/KOKKOS/fix_wall_region_kokkos.cpp index cc8d876d043..167ac394660 100644 --- a/src/KOKKOS/fix_wall_region_kokkos.cpp +++ b/src/KOKKOS/fix_wall_region_kokkos.cpp @@ -28,6 +28,8 @@ #include "region_block_kokkos.h" #include "region_sphere_kokkos.h" +#include + using namespace LAMMPS_NS; using namespace MathSpecialKokkos; @@ -78,9 +80,12 @@ void FixWallRegionKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); @@ -96,10 +101,6 @@ void FixWallRegionKokkos::post_force(int vflag) region->prematch(); - // virial setup - - v_init(vflag); - // region->match() ensures particle is in region or on surface, else error // if returned contact dist r = 0, is on surface, also an error // in COLLOID case, r <= radius is an error @@ -164,14 +165,21 @@ void FixWallRegionKokkos::wall_particle(T regionKK, const int i, val else tooclose = 0.0; - int n = regionKK->surface_kokkos(static_cast(d_x(i,0)), static_cast(d_x(i,1)), static_cast(d_x(i,2)), cutoff); + // the contact list lives on the stack, so that concurrently running + // threads do not overwrite each other's contacts + + Region::Contact contact[std::remove_pointer_t::MAXCONTACT]; + + int n = regionKK->surface_kokkos(static_cast(d_x(i,0)), + static_cast(d_x(i,1)), + static_cast(d_x(i,2)), cutoff, contact); for ( int m = 0; m < n; m++) { - KK_FLOAT r = static_cast(regionKK->d_contact[m].r); - KK_FLOAT delx = static_cast(regionKK->d_contact[m].delx); - KK_FLOAT dely = static_cast(regionKK->d_contact[m].dely); - KK_FLOAT delz = static_cast(regionKK->d_contact[m].delz); + KK_FLOAT r = static_cast(contact[m].r); + KK_FLOAT delx = static_cast(contact[m].delx); + KK_FLOAT dely = static_cast(contact[m].dely); + KK_FLOAT delz = static_cast(contact[m].delz); if (r <= tooclose) Kokkos::abort("Particle outside surface of region used in fix wall/region"); diff --git a/src/KOKKOS/grid3d_kokkos.cpp b/src/KOKKOS/grid3d_kokkos.cpp index b0b20020944..e200a1e8b80 100644 --- a/src/KOKKOS/grid3d_kokkos.cpp +++ b/src/KOKKOS/grid3d_kokkos.cpp @@ -711,7 +711,7 @@ void Grid3dKokkos:: forward_comm_kspace_tiled(KSpace *kspace, int which, int nper, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { - int i,m,offset; + int m,offset; KokkosBaseFFT* kspaceKKBase = dynamic_cast(kspace); FFT_SCALAR* buf1; @@ -755,14 +755,19 @@ forward_comm_kspace_tiled(KSpace *kspace, int which, int nper, // unpack all received data - for (i = 0; i < nrecv; i++) { - MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE); + // every receive writes into its own stretch of the one host buffer, and the + // copy to the device takes the whole buffer, so it cannot be started while + // any of them is still in flight: it would read the stretches that have not + // arrived and copy them over the device side. Collect them all first, copy + // once, then unpack. The brick path already waits before it copies. + MPI_Waitall(nrecv,requests,MPI_STATUSES_IGNORE); - if (!lmp->kokkos->gpu_aware_flag) { - k_buf2.modify_host(); - k_buf2.sync(); - } + if (!lmp->kokkos->gpu_aware_flag) { + k_buf2.modify_host(); + k_buf2.sync(); + } + for (m = 0; m < nrecv; m++) { offset = nper * recv[m].offset; kspaceKKBase->unpack_forward_grid_kokkos(which,k_buf2,offset, recv[m].nunpack,k_recv_unpacklist,m); @@ -852,7 +857,7 @@ void Grid3dKokkos:: reverse_comm_kspace_tiled(KSpace *kspace, int which, int nper, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { - int i,m,offset; + int m,offset; KokkosBaseFFT* kspaceKKBase = dynamic_cast(kspace); @@ -896,14 +901,16 @@ reverse_comm_kspace_tiled(KSpace *kspace, int which, int nper, } // unpack all received data - for (i = 0; i < nsend; i++) { - MPI_Waitany(nsend,requests,&m,MPI_STATUS_IGNORE); + // as in the forward direction: one buffer, one copy, and only once every + // receive into it has landed + MPI_Waitall(nsend,requests,MPI_STATUSES_IGNORE); - if (!lmp->kokkos->gpu_aware_flag) { - k_buf2.modify_host(); - k_buf2.sync(); - } + if (!lmp->kokkos->gpu_aware_flag) { + k_buf2.modify_host(); + k_buf2.sync(); + } + for (m = 0; m < nsend; m++) { offset = nper * send[m].offset; kspaceKKBase->unpack_reverse_grid_kokkos(which,k_buf2,offset, send[m].npack,k_send_packlist,m); diff --git a/src/KOKKOS/improper_class2_kokkos.cpp b/src/KOKKOS/improper_class2_kokkos.cpp index c41e57b8ccf..ca34eba536e 100644 --- a/src/KOKKOS/improper_class2_kokkos.cpp +++ b/src/KOKKOS/improper_class2_kokkos.cpp @@ -82,7 +82,6 @@ void ImproperClass2Kokkos::compute(int eflag_in, int vflag_in) d_vatom = k_vatom.template view(); } - atomKK->sync(execution_space,datamask_read); k_k0.template sync(); k_chi0.template sync(); k_aa_k1.template sync(); @@ -95,8 +94,6 @@ void ImproperClass2Kokkos::compute(int eflag_in, int vflag_in) k_setflag_i.template sync(); k_setflag_aa.template sync(); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_cossq_kokkos.cpp b/src/KOKKOS/improper_cossq_kokkos.cpp index f633eb69d33..be7f1a7fc70 100644 --- a/src/KOKKOS/improper_cossq_kokkos.cpp +++ b/src/KOKKOS/improper_cossq_kokkos.cpp @@ -88,9 +88,6 @@ void ImproperCossqKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_cvff_kokkos.cpp b/src/KOKKOS/improper_cvff_kokkos.cpp index e734a6d6855..c7020cb0f9e 100644 --- a/src/KOKKOS/improper_cvff_kokkos.cpp +++ b/src/KOKKOS/improper_cvff_kokkos.cpp @@ -88,12 +88,9 @@ void ImproperCvffKokkos::compute(int eflag_in, int vflag_in) } else Kokkos::deep_copy(d_vatom,0.0); } - atomKK->sync(execution_space,datamask_read); k_k.template sync(); k_sign.template sync(); k_multiplicity.template sync(); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_distance_kokkos.cpp b/src/KOKKOS/improper_distance_kokkos.cpp index 41ce2c82c4b..68a00c0c02f 100644 --- a/src/KOKKOS/improper_distance_kokkos.cpp +++ b/src/KOKKOS/improper_distance_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperDistanceKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_distharm_kokkos.cpp b/src/KOKKOS/improper_distharm_kokkos.cpp index 8ae6cddbcb0..f837f398433 100644 --- a/src/KOKKOS/improper_distharm_kokkos.cpp +++ b/src/KOKKOS/improper_distharm_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperDistHarmKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_fourier_kokkos.cpp b/src/KOKKOS/improper_fourier_kokkos.cpp index 367669321e3..689bf3c84ab 100644 --- a/src/KOKKOS/improper_fourier_kokkos.cpp +++ b/src/KOKKOS/improper_fourier_kokkos.cpp @@ -87,9 +87,6 @@ void ImproperFourierKokkos::compute(int eflag_in, int vflag_in) k_C1.template sync(); k_C2.template sync(); k_all.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_harmonic_kokkos.cpp b/src/KOKKOS/improper_harmonic_kokkos.cpp index b92a2afbb88..6f7ba766925 100644 --- a/src/KOKKOS/improper_harmonic_kokkos.cpp +++ b/src/KOKKOS/improper_harmonic_kokkos.cpp @@ -88,11 +88,8 @@ void ImproperHarmonicKokkos::compute(int eflag_in, int vflag_in) } else Kokkos::deep_copy(d_vatom,0.0); } - atomKK->sync(execution_space,datamask_read); k_k.template sync(); k_chi.template sync(); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp b/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp index a51ed5eea1d..47901f2d947 100644 --- a/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp +++ b/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperInversionHarmonicKokkos::compute(int eflag_in, int vfla k_kw.template sync(); k_w0.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_ring_kokkos.cpp b/src/KOKKOS/improper_ring_kokkos.cpp index 58071b7fdef..82b0b416851 100644 --- a/src/KOKKOS/improper_ring_kokkos.cpp +++ b/src/KOKKOS/improper_ring_kokkos.cpp @@ -83,9 +83,6 @@ void ImproperRingKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_sqdistharm_kokkos.cpp b/src/KOKKOS/improper_sqdistharm_kokkos.cpp index 87dccfb8393..b238e2e4d5e 100644 --- a/src/KOKKOS/improper_sqdistharm_kokkos.cpp +++ b/src/KOKKOS/improper_sqdistharm_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperSQDistHarmKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_umbrella_kokkos.cpp b/src/KOKKOS/improper_umbrella_kokkos.cpp index def6da5ed0c..003fdbe587b 100644 --- a/src/KOKKOS/improper_umbrella_kokkos.cpp +++ b/src/KOKKOS/improper_umbrella_kokkos.cpp @@ -85,9 +85,6 @@ void ImproperUmbrellaKokkos::compute(int eflag_in, int vflag_in) k_kw.template sync(); k_w0.template sync(); k_C.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/min_fire_kokkos.cpp b/src/KOKKOS/min_fire_kokkos.cpp index 84f95aac1e0..944db143e10 100644 --- a/src/KOKKOS/min_fire_kokkos.cpp +++ b/src/KOKKOS/min_fire_kokkos.cpp @@ -120,9 +120,28 @@ int MinFireKokkos::run_iterate(int maxiter) { auto l_type = atomKK->k_type.view_device(); int nlocal = atom->nlocal; + // energy_force() runs the communication, the neighbor build and the fixes. + // Any of those can leave the newest per-atom data on the host -- a fix + // without KOKKOS support writes there -- and can grow or reorder the arrays, + // so the data and the views both have to be taken again after every call. + // Without this the kernels below integrate a device copy that the host has + // since overtaken, and the modified() that follows finds both sides claimed. + + auto refresh = [&]() { + atomKK->sync(Device, X_MASK | V_MASK | F_MASK | RMASS_MASK | TYPE_MASK); + l_x = atomKK->k_x.view_device(); + l_v = atomKK->k_v.view_device(); + l_f = atomKK->k_f.view_device(); + l_rmass = atomKK->k_rmass.view_device(); + l_mass = atomKK->k_mass.view_device(); + l_type = atomKK->k_type.view_device(); + nlocal = atom->nlocal; + }; + if constexpr (INTEGRATOR == LEAPFROG) { energy_force(0); neval++; + refresh(); double dtf = -0.5 * dt * force->ftm2v; const KK_FLOAT dtf_kk = static_cast(dtf); Kokkos::parallel_for("min_fire/leapfrog_init", atom->nlocal, LAMMPS_LAMBDA(const int i) { @@ -131,11 +150,14 @@ int MinFireKokkos::run_iterate(int maxiter) { l_v(i,1) = dtfm * static_cast(l_f(i,1)); l_v(i,2) = dtfm * static_cast(l_f(i,2)); }); + atomKK->modified(Device, V_MASK); } for (int iter = 0; iter < maxiter; iter++) { if (timer->check_timeout(niter)) return TIMEOUT; + refresh(); + bigint ntimestep = ++update->ntimestep; niter++; @@ -210,12 +232,14 @@ int MinFireKokkos::run_iterate(int maxiter) { } l_v(i,0) = l_v(i,1) = l_v(i,2) = 0.0; }); + atomKK->modified(Device, X_MASK | V_MASK); flagv0 = 1; } if (!ABCFLAG && flagv0) { energy_force(0); neval++; + refresh(); double dtf_init = dt * force->ftm2v; const KK_FLOAT dtf_init_kk = static_cast(dtf_init); Kokkos::parallel_for("min_fire/v_init", nlocal, LAMMPS_LAMBDA(const int i) { @@ -224,6 +248,7 @@ int MinFireKokkos::run_iterate(int maxiter) { l_v(i,1) = dtfm * static_cast(l_f(i,1)); l_v(i,2) = dtfm * static_cast(l_f(i,2)); }); + atomKK->modified(Device, V_MASK); } // cannot use "if constexpr" below because CUDA device lambdas @@ -249,6 +274,7 @@ int MinFireKokkos::run_iterate(int maxiter) { Kokkos::parallel_for("min_fire/final_v_zero", nlocal, LAMMPS_LAMBDA(const int i) { l_v(i,0) = l_v(i,1) = l_v(i,2) = 0.0; }); + atomKK->modified(Device, V_MASK); } KK_FLOAT dtf_final = dtv * static_cast(force->ftm2v); @@ -318,9 +344,9 @@ int MinFireKokkos::run_iterate(int maxiter) { eprevious = ecurrent; ecurrent = energy_force(0); neval++; + refresh(); if constexpr (INTEGRATOR == VERLET) { - atomKK->sync(Device, V_MASK | F_MASK); Kokkos::parallel_for("min_fire/verlet_v_final", nlocal, LAMMPS_LAMBDA(const int i) { KK_FLOAT dtfm_half = dtf_half / (l_rmass.data() ? l_rmass(i) : l_mass(l_type(i))); l_v(i,0) += dtfm_half * static_cast(l_f(i,0)); diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 26ee88ff513..4ff0fc1b48a 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -53,7 +53,18 @@ void ModifyKokkos::setup(int vflag) } } - for (int i = 0; i < ncompute; i++) compute[i]->setup(); + // the computes get the same treatment as the fixes above: several of them + // read the per-atom arrays in setup(), through the plain pointers when they + // have no KOKKOS version, and were reaching them without a transfer + + for (int i = 0; i < ncompute; i++) { + atomKK->sync(compute[i]->execution_space,compute[i]->datamask_read); + int prev_auto_sync = lmp->kokkos->auto_sync; + if (!compute[i]->kokkosable) lmp->kokkos->auto_sync = 1; + compute[i]->setup(); + lmp->kokkos->auto_sync = prev_auto_sync; + atomKK->modified(compute[i]->execution_space,compute[i]->datamask_modify); + } if (update->whichflag == 1) for (int i = 0; i < nfix; i++) { @@ -451,6 +462,8 @@ double ModifyKokkos::energy_couple() { double energy = 0.0; for (int i = 0; i < n_energy_couple; i++) { + atomKK->sync(fix[list_energy_couple[i]]->execution_space, + fix[list_energy_couple[i]]->datamask_read); int prev_auto_sync = lmp->kokkos->auto_sync; if (!fix[list_energy_couple[i]]->kokkosable) lmp->kokkos->auto_sync = 1; energy += fix[list_energy_couple[i]]->compute_scalar(); @@ -471,6 +484,8 @@ double ModifyKokkos::energy_global() { double energy = 0.0; for (int i = 0; i < n_energy_global; i++) { + atomKK->sync(fix[list_energy_global[i]]->execution_space, + fix[list_energy_global[i]]->datamask_read); int prev_auto_sync = lmp->kokkos->auto_sync; if (!fix[list_energy_global[i]]->kokkosable) lmp->kokkos->auto_sync = 1; energy += fix[list_energy_global[i]]->compute_scalar(); diff --git a/src/KOKKOS/nbin_kokkos.cpp b/src/KOKKOS/nbin_kokkos.cpp index 694f4334382..3652b0bf338 100644 --- a/src/KOKKOS/nbin_kokkos.cpp +++ b/src/KOKKOS/nbin_kokkos.cpp @@ -122,7 +122,32 @@ void NBinKokkos::bin_atoms() Kokkos::deep_copy(h_resize, d_resize); if (h_resize()) { - atoms_per_bin += 16; + // A bin overflowed its capacity. bincount now holds the true + // occupancy of every bin, because the atomic increment in + // binatomsItem() runs for every binned atom whether or not the bin + // was already full. Size atoms_per_bin from the actual maximum + // occupancy in a single step instead of growing by a fixed increment + // and re-binning all atoms once per increment. The latter costs + // O(max bin occupancy) reallocations and re-bins and dominates + // neighbor setup for skewed distributions, e.g. the large cutoff bins + // used by the KOKKOS package on GPUs. + + auto d_bincount = k_bincount.view(); + int max_bincount = 0; + Kokkos::parallel_reduce(mbins, + LAMMPS_LAMBDA(const int i, int &max_val) { + max_val = MAX(max_val,d_bincount[i]); + },Kokkos::Max(max_bincount)); + + // grow to the true maximum occupancy plus ~10% headroom (at least 16) + // so small density fluctuations on later steps do not immediately + // force another regrow. Reaching this branch means a bin overflowed, + // so max_bincount > atoms_per_bin and the new capacity strictly + // exceeds both the old one and the true occupancy: the next pass + // cannot overflow, bounding the loop at one more re-bin. + + atoms_per_bin = max_bincount + MAX(16,max_bincount/10); + k_bins = DAT::tdual_int_2d("Neighbor::bins", mbins, atoms_per_bin); bins = k_bins.view(); c_bins = bins; @@ -167,6 +192,7 @@ void NBinKokkos::binatomsItem(const int &i) const Kokkos::abort("Atom outside of neighbor bin range - simulation unstable"); atom2bin(i) = ibin; + const int ac = Kokkos::atomic_fetch_add(&bincount[ibin], (int)1); if (ac < (int)bins.extent(1)) { bins(ibin, ac) = i; diff --git a/src/KOKKOS/neigh_bond_kokkos.cpp b/src/KOKKOS/neigh_bond_kokkos.cpp index b1205895eae..f083afdcb56 100644 --- a/src/KOKKOS/neigh_bond_kokkos.cpp +++ b/src/KOKKOS/neigh_bond_kokkos.cpp @@ -66,6 +66,8 @@ NeighBondKokkos::NeighBondKokkos(LAMMPS *lmp) : Pointers(lmp) maxangle = 0; maxdihedral = 0; maximproper = 0; + + copymode = 0; } /* ---------------------------------------------------------------------- */ @@ -73,6 +75,12 @@ NeighBondKokkos::NeighBondKokkos(LAMMPS *lmp) : Pointers(lmp) template NeighBondKokkos::~NeighBondKokkos() { + // Kokkos hands the loop bodies below a copy of this class as the functor and + // destroys that copy when the loop is done. The topology lists belong to + // Neighbor, not to the copy, so only the original may release them. + + if (copymode) return; + memoryKK->destroy_kokkos(k_bondlist,neighbor->bondlist); memoryKK->destroy_kokkos(k_anglelist,neighbor->anglelist); memoryKK->destroy_kokkos(k_dihedrallist,neighbor->dihedrallist); @@ -124,6 +132,12 @@ void NeighBondKokkos::init_topology_kk() { int i,m; int bond_off = 0; int angle_off = 0; + // keep this list the same as the one in Neighbor::init_topology(): a fix + // that turns bonds off by making their type negative has to be named here, + // because it does so after this decision is taken and the scan below cannot + // see it yet. Without the name the all variant is chosen, and that one + // copies the type into the list without looking at its sign. + for (const auto &ifix : modify->get_fix_list()) if (utils::strmatch(ifix->style,"^shake") || utils::strmatch(ifix->style,"^rattle") || utils::strmatch(ifix->style,"^ilves")) @@ -241,6 +255,11 @@ template void NeighBondKokkos::bond_all() { atomKK->sync(execution_space, BOND_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_bondlist.clear_sync_state(); v_bondlist = k_bondlist.view(); num_bond = atomKK->k_num_bond.view(); bond_atom = atomKK->k_bond_atom.view(); @@ -256,7 +275,9 @@ void NeighBondKokkos::bond_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -322,6 +343,11 @@ template void NeighBondKokkos::bond_partial() { atomKK->sync(execution_space, BOND_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_bondlist.clear_sync_state(); v_bondlist = k_bondlist.view(); num_bond = atomKK->k_num_bond.view(); bond_atom = atomKK->k_bond_atom.view(); @@ -337,7 +363,9 @@ void NeighBondKokkos::bond_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -399,7 +427,9 @@ void NeighBondKokkos::bond_check() atomKK->sync(execution_space, X_MASK); k_bondlist.sync(); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,neighbor->nbondlist),*this,flag); + copymode = 0; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); @@ -427,6 +457,11 @@ template void NeighBondKokkos::angle_all() { atomKK->sync(execution_space, ANGLE_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_anglelist.clear_sync_state(); v_anglelist = k_anglelist.view(); num_angle = atomKK->k_num_angle.view(); angle_atom1 = atomKK->k_angle_atom1.view(); @@ -444,7 +479,9 @@ void NeighBondKokkos::angle_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -514,6 +551,11 @@ template void NeighBondKokkos::angle_partial() { atomKK->sync(execution_space, ANGLE_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_anglelist.clear_sync_state(); v_anglelist = k_anglelist.view(); num_angle = atomKK->k_num_angle.view(); angle_atom1 = atomKK->k_angle_atom1.view(); @@ -531,7 +573,9 @@ void NeighBondKokkos::angle_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -601,7 +645,9 @@ void NeighBondKokkos::angle_check() atomKK->sync(execution_space, X_MASK); k_anglelist.sync(); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,neighbor->nanglelist),*this,flag); + copymode = 0; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); @@ -640,6 +686,11 @@ template void NeighBondKokkos::dihedral_all() { atomKK->sync(execution_space, DIHEDRAL_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_dihedrallist.clear_sync_state(); v_dihedrallist = k_dihedrallist.view(); num_dihedral = atomKK->k_num_dihedral.view(); dihedral_atom1 = atomKK->k_dihedral_atom1.view(); @@ -658,7 +709,9 @@ void NeighBondKokkos::dihedral_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -732,6 +785,11 @@ template void NeighBondKokkos::dihedral_partial() { atomKK->sync(execution_space, DIHEDRAL_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_dihedrallist.clear_sync_state(); v_dihedrallist = k_dihedrallist.view(); num_dihedral = atomKK->k_num_dihedral.view(); dihedral_atom1 = atomKK->k_dihedral_atom1.view(); @@ -750,7 +808,9 @@ void NeighBondKokkos::dihedral_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -825,7 +885,9 @@ void NeighBondKokkos::dihedral_check(int nlist, typename AT::t_int_2 atomKK->sync(execution_space, X_MASK); k_dihedrallist.sync(); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlist),*this,flag); + copymode = 0; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); @@ -881,6 +943,11 @@ template void NeighBondKokkos::improper_all() { atomKK->sync(execution_space, IMPROPER_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_improperlist.clear_sync_state(); v_improperlist = k_improperlist.view(); num_improper = atomKK->k_num_improper.view(); improper_atom1 = atomKK->k_improper_atom1.view(); @@ -899,7 +966,9 @@ void NeighBondKokkos::improper_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -973,6 +1042,11 @@ template void NeighBondKokkos::improper_partial() { atomKK->sync(execution_space, IMPROPER_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_improperlist.clear_sync_state(); v_improperlist = k_improperlist.view(); num_improper = atomKK->k_num_improper.view(); improper_atom1 = atomKK->k_improper_atom1.view(); @@ -991,7 +1065,9 @@ void NeighBondKokkos::improper_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); diff --git a/src/KOKKOS/neigh_bond_kokkos.h b/src/KOKKOS/neigh_bond_kokkos.h index 38208b99f64..daea6b92db4 100644 --- a/src/KOKKOS/neigh_bond_kokkos.h +++ b/src/KOKKOS/neigh_bond_kokkos.h @@ -90,6 +90,7 @@ class NeighBondKokkos : protected Pointers { int maxbond,maxangle,maxdihedral,maximproper; // size of bond lists int me,nprocs; + int copymode; // 1 while a Kokkos loop holds a copy private: int map_style; diff --git a/src/KOKKOS/neighbor_kokkos.cpp b/src/KOKKOS/neighbor_kokkos.cpp index 57db7b54b9a..dd7e9922d6c 100644 --- a/src/KOKKOS/neighbor_kokkos.cpp +++ b/src/KOKKOS/neighbor_kokkos.cpp @@ -420,17 +420,6 @@ void NeighborKokkos::build_topology() { k_dihedrallist = neighbond_device.k_dihedrallist; k_improperlist = neighbond_device.k_improperlist; - // Transfer topology neighbor lists to Host for non-Kokkos styles - - if (force->bond && force->bond->execution_space == Host) - k_bondlist.sync_host(); - if (force->angle && force->angle->execution_space == Host) - k_anglelist.sync_host(); - if (force->dihedral && force->dihedral->execution_space == Host) - k_dihedrallist.sync_host(); - if (force->improper && force->improper->execution_space == Host) - k_improperlist.sync_host(); - } else { neighbond_host.build_topology_kk(); @@ -439,4 +428,16 @@ void NeighborKokkos::build_topology() { k_dihedrallist = neighbond_host.k_dihedrallist; k_improperlist = neighbond_host.k_improperlist; } + + // transfer topology neighbor lists to the host for non-Kokkos styles, + // which read them through the plain pointers in Neighbor + + if (force->bond && force->bond->execution_space == Host) + k_bondlist.sync_host(); + if (force->angle && force->angle->execution_space == Host) + k_anglelist.sync_host(); + if (force->dihedral && force->dihedral->execution_space == Host) + k_dihedrallist.sync_host(); + if (force->improper && force->improper->execution_space == Host) + k_improperlist.sync_host(); } diff --git a/src/KOKKOS/pair_buck_coul_long_kokkos.cpp b/src/KOKKOS/pair_buck_coul_long_kokkos.cpp index 4ee2bb0526b..927335c9eea 100644 --- a/src/KOKKOS/pair_buck_coul_long_kokkos.cpp +++ b/src/KOKKOS/pair_buck_coul_long_kokkos.cpp @@ -142,6 +142,16 @@ void PairBuckCoulLongKokkos::compute(int eflag_in, int vflag_in) virial[5] += static_cast(ev.v[5]); } + if (eflag_atom) { + k_eatom.template modify(); + k_eatom.sync_host(); + } + + if (vflag_atom) { + k_vatom.template modify(); + k_vatom.sync_host(); + } + if (vflag_fdotr) pair_virial_fdotr_compute(this); copymode = 0; diff --git a/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp b/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp index e0291f72476..f1f50fbf150 100644 --- a/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp +++ b/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp @@ -164,8 +164,6 @@ void PairDPDfdtEnergyKokkos::compute(int eflag_in, int vflag_in) k_cutsq.template sync(); k_params.template sync(); atomKK->sync(execution_space,X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK); - if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); - else atomKK->modified(execution_space,F_MASK); special_lj[0] = static_cast(force->special_lj[0]); special_lj[1] = static_cast(force->special_lj[1]); @@ -328,8 +326,26 @@ void PairDPDfdtEnergyKokkos::compute(int eflag_in, int vflag_in) k_duMech.template modify(); k_duMech.sync_host(); comm->reverse_comm(this); + + // the reverse communication adds the ghost contributions through the plain + // host arrays, so claim that write: without it the device copy keeps the + // pre-communication values and the sync in fix dpd/energy/kk has nothing + // left to copy + + k_duCond.modify_host(); + k_duMech.modify_host(); } + + // claim the forces here rather than before the kernels above: this style + // declares an empty datamask, so this is the only claim they get, and the + // reverse communication of the energy changes syncs the host in between, + // which would take a claim made up there with it and leave the forces the + // kernels wrote unclaimed + + if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); + else atomKK->modified(execution_space,F_MASK); + if (eflag_global) eng_vdwl += static_cast(ev.evdwl); if (vflag_global) { virial[0] += static_cast(ev.v[0]); diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index 21d5d786e46..9f707c25d54 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -28,6 +28,7 @@ #include "memory_kokkos.h" #include "neighbor_kokkos.h" #include "neigh_request.h" +#include "utils.h" #include "ace-evaluator/ace_version.h" #include "ace-evaluator/ace_radial.h" @@ -64,6 +65,10 @@ PairPACEExtrapolationKokkos::PairPACEExtrapolationKokkos(LAMMPS *lmp datamask_modify = EMPTY_MASK; host_flag = (execution_space == HostKK); + + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + neigh_scratch_level = 0; + neigh_scratch_warned = 0; } /* ---------------------------------------------------------------------- @@ -501,6 +506,76 @@ double PairPACEExtrapolationKokkos::init_one(int i, int j) return cutone; } +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +template +void PairPACEExtrapolationKokkos::settings(int narg, char **arg) +{ + // intercept the KOKKOS-only "neigh" keyword, which selects the team scratch + // memory level used to build the short neighbor list, then forward the + // remaining keywords to the CPU base class + + auto base_arg = new char*[narg]; + int base_narg = 0; + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg], "neigh") == 0) { + if (iarg+2 > narg) + utils::missing_cmd_args(FLERR, "pair_style pace/extrapolation neigh", error); + if (strcmp(arg[iarg+1], "auto") == 0) + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + else if (strcmp(arg[iarg+1], "shared") == 0) + neigh_scratch_request = NEIGH_SCRATCH_SHARED; + else if (strcmp(arg[iarg+1], "global") == 0) + neigh_scratch_request = NEIGH_SCRATCH_GLOBAL; + else + error->all(FLERR, "Unknown pair_style pace/extrapolation neigh keyword: {}", arg[iarg+1]); + iarg += 2; + } else { + base_arg[base_narg++] = arg[iarg]; + iarg++; + } + } + + PairPACEExtrapolation::settings(base_narg, base_arg); + + delete[] base_arg; +} + +/* ---------------------------------------------------------------------- + select the team scratch memory level for the ComputeNeigh short neighbor + list build; falls back from level 0 (fast on-chip shared memory) to level 1 + (global memory) when the request does not fit into the available shared + memory, unless the user forced a level +------------------------------------------------------------------------- */ + +template +int PairPACEExtrapolationKokkos::neigh_scratch_level_select(int scratch_size, int max_level0) +{ + // honor an explicit user request + if (neigh_scratch_request == NEIGH_SCRATCH_SHARED) return 0; + if (neigh_scratch_request == NEIGH_SCRATCH_GLOBAL) return 1; + + // automatic: use fast level-0 (shared) scratch when it fits, otherwise fall + // back to level-1 (global) scratch. max_level0 is queried from Kokkos rather + // than hard-coded, so larger shared-memory limits (e.g. the opt-in >48 KiB + // shared memory available in newer Kokkos) are used automatically. + if (scratch_size <= max_level0) return 0; + + if (!neigh_scratch_warned && comm->me == 0) { + error->warning(FLERR, + "Pair pace/extrapolation/kk short neighbor list needs {} bytes of team " + "scratch memory but only {} bytes of on-chip (level-0) shared memory are " + "available; falling back to slower global (level-1) memory. Reduce the " + "neighbor count or use the pair_style 'neigh global' keyword to silence " + "this warning.", scratch_size, max_level0); + neigh_scratch_warned = 1; + } + return 1; +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ @@ -676,7 +751,16 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in check_team_size_for(chunk_size,team_size,vector_length); int scratch_size = scratch_size_helper(team_size * maxneigh); typename Kokkos::TeamPolicy policy_neigh(chunk_size,team_size,vector_length); - policy_neigh = policy_neigh.set_scratch_size(0, Kokkos::PerTeam(scratch_size)); + + // The ComputeNeigh kernel caches the short neighbor list in team scratch + // memory. On GPUs level-0 scratch is fast on-chip shared memory but is a + // scarce resource: with many neighbors and/or atom types the request can + // exceed what the device provides and abort the run (see + // https://github.com/lammps/lammps/issues/5063). Query the level-0 limit + // from Kokkos (never hard-coded) and transparently fall back to level-1 + // (global memory) scratch when the request does not fit. + neigh_scratch_level = neigh_scratch_level_select(scratch_size, policy_neigh.scratch_size_max(0)); + policy_neigh = policy_neigh.set_scratch_size(neigh_scratch_level, Kokkos::PerTeam(scratch_size)); Kokkos::parallel_for("ComputeNeigh",policy_neigh,*this); } @@ -835,7 +919,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig // If it is, inside is assigned to 1, otherwise -1 const int team_rank = team.team_rank(); const int scratch_shift = team_rank * maxneigh; // offset into pointer for entire team - int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), 0) + scratch_shift; + int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), neigh_scratch_level) + scratch_shift; // loop over list of all neighbors within force cutoff // distsq[] = distance sq to each diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.h b/src/KOKKOS/pair_pace_extrapolation_kokkos.h index 4aec48a7c3f..bc76c7fc0bb 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.h +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.h @@ -56,6 +56,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { ~PairPACEExtrapolationKokkos() override; void compute(int, int) override; + void settings(int, char **) override; void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; @@ -110,6 +111,19 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max, idx_sph_max; int host_flag; + // team scratch memory level used by the ComputeNeigh short neighbor list build: + // NEIGH_SCRATCH_AUTO - automatically use level 0 (fast on-chip shared + // memory) when it fits, else fall back to level 1 + // (global memory) + // NEIGH_SCRATCH_SHARED - always use level 0 + // NEIGH_SCRATCH_GLOBAL - always use level 1 + enum { NEIGH_SCRATCH_AUTO = 0, NEIGH_SCRATCH_SHARED, NEIGH_SCRATCH_GLOBAL }; + int neigh_scratch_request; // user preference (pair_style "neigh" keyword) + int neigh_scratch_level; // level actually used by ComputeNeigh (0 or 1) + int neigh_scratch_warned; // whether the auto-fallback warning was printed + + int neigh_scratch_level_select(int scratch_size, int max_level0); + int eflag, vflag; int neighflag, max_ndensity; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index f954fd4b69c..e72af91da85 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -28,6 +28,7 @@ #include "memory_kokkos.h" #include "neighbor_kokkos.h" #include "neigh_request.h" +#include "utils.h" #include "ace-evaluator/ace_version.h" #include "ace-evaluator/ace_radial.h" @@ -65,6 +66,10 @@ PairPACEKokkos::PairPACEKokkos(LAMMPS *lmp) : PairPACE(lmp) datamask_modify = EMPTY_MASK; host_flag = (execution_space == HostKK); + + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + neigh_scratch_level = 0; + neigh_scratch_warned = 0; } /* ---------------------------------------------------------------------- @@ -473,6 +478,76 @@ double PairPACEKokkos::init_one(int i, int j) return cutone; } +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +template +void PairPACEKokkos::settings(int narg, char **arg) +{ + // intercept the KOKKOS-only "neigh" keyword, which selects the team scratch + // memory level used to build the short neighbor list, then forward the + // remaining keywords to the CPU base class + + auto base_arg = new char*[narg]; + int base_narg = 0; + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg], "neigh") == 0) { + if (iarg+2 > narg) + utils::missing_cmd_args(FLERR, "pair_style pace neigh", error); + if (strcmp(arg[iarg+1], "auto") == 0) + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + else if (strcmp(arg[iarg+1], "shared") == 0) + neigh_scratch_request = NEIGH_SCRATCH_SHARED; + else if (strcmp(arg[iarg+1], "global") == 0) + neigh_scratch_request = NEIGH_SCRATCH_GLOBAL; + else + error->all(FLERR, "Unknown pair_style pace neigh keyword: {}", arg[iarg+1]); + iarg += 2; + } else { + base_arg[base_narg++] = arg[iarg]; + iarg++; + } + } + + PairPACE::settings(base_narg, base_arg); + + delete[] base_arg; +} + +/* ---------------------------------------------------------------------- + select the team scratch memory level for the ComputeNeigh short neighbor + list build; falls back from level 0 (fast on-chip shared memory) to level 1 + (global memory) when the request does not fit into the available shared + memory, unless the user forced a level +------------------------------------------------------------------------- */ + +template +int PairPACEKokkos::neigh_scratch_level_select(int scratch_size, int max_level0) +{ + // honor an explicit user request + if (neigh_scratch_request == NEIGH_SCRATCH_SHARED) return 0; + if (neigh_scratch_request == NEIGH_SCRATCH_GLOBAL) return 1; + + // automatic: use fast level-0 (shared) scratch when it fits, otherwise fall + // back to level-1 (global) scratch. max_level0 is queried from Kokkos rather + // than hard-coded, so larger shared-memory limits (e.g. the opt-in >48 KiB + // shared memory available in newer Kokkos) are used automatically. + if (scratch_size <= max_level0) return 0; + + if (!neigh_scratch_warned && comm->me == 0) { + error->warning(FLERR, + "Pair pace/kk short neighbor list needs {} bytes of team scratch memory " + "but only {} bytes of on-chip (level-0) shared memory are available; " + "falling back to slower global (level-1) memory. Reduce the neighbor " + "count or use the pair_style 'neigh global' keyword to silence this " + "warning.", scratch_size, max_level0); + neigh_scratch_warned = 1; + } + return 1; +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ @@ -636,7 +711,16 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) check_team_size_for(chunk_size,team_size,vector_length); int scratch_size = scratch_size_helper(team_size * maxneigh); typename Kokkos::TeamPolicy policy_neigh(chunk_size,team_size,vector_length); - policy_neigh = policy_neigh.set_scratch_size(0, Kokkos::PerTeam(scratch_size)); + + // The ComputeNeigh kernel caches the short neighbor list in team scratch + // memory. On GPUs level-0 scratch is fast on-chip shared memory but is a + // scarce resource: with many neighbors and/or atom types the request can + // exceed what the device provides and abort the run (see + // https://github.com/lammps/lammps/issues/5063). Query the level-0 limit + // from Kokkos (never hard-coded) and transparently fall back to level-1 + // (global memory) scratch when the request does not fit. + neigh_scratch_level = neigh_scratch_level_select(scratch_size, policy_neigh.scratch_size_max(0)); + policy_neigh = policy_neigh.set_scratch_size(neigh_scratch_level, Kokkos::PerTeam(scratch_size)); Kokkos::parallel_for("ComputeNeigh",policy_neigh,*this); } @@ -781,7 +865,7 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen // If it is, inside is assigned to 1, otherwise -1 const int team_rank = team.team_rank(); const int scratch_shift = team_rank * maxneigh; // offset into pointer for entire team - int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), 0) + scratch_shift; + int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), neigh_scratch_level) + scratch_shift; // loop over list of all neighbors within force cutoff // distsq[] = distance sq to each diff --git a/src/KOKKOS/pair_pace_kokkos.h b/src/KOKKOS/pair_pace_kokkos.h index cd3cd6b037f..f1203b8386c 100644 --- a/src/KOKKOS/pair_pace_kokkos.h +++ b/src/KOKKOS/pair_pace_kokkos.h @@ -55,6 +55,7 @@ class PairPACEKokkos : public PairPACE { ~PairPACEKokkos() override; void compute(int, int) override; + void settings(int, char **) override; void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; @@ -105,6 +106,19 @@ class PairPACEKokkos : public PairPACE { int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, idx_sph_max; int host_flag; + // team scratch memory level used by the ComputeNeigh short neighbor list build: + // NEIGH_SCRATCH_AUTO - automatically use level 0 (fast on-chip shared + // memory) when it fits, else fall back to level 1 + // (global memory) + // NEIGH_SCRATCH_SHARED - always use level 0 + // NEIGH_SCRATCH_GLOBAL - always use level 1 + enum { NEIGH_SCRATCH_AUTO = 0, NEIGH_SCRATCH_SHARED, NEIGH_SCRATCH_GLOBAL }; + int neigh_scratch_request; // user preference (pair_style "neigh" keyword) + int neigh_scratch_level; // level actually used by ComputeNeigh (0 or 1) + int neigh_scratch_warned; // whether the auto-fallback warning was printed + + int neigh_scratch_level_select(int scratch_size, int max_level0); + int eflag, vflag; int neighflag, max_ndensity; diff --git a/src/KOKKOS/pair_tip4p_kokkos.h b/src/KOKKOS/pair_tip4p_kokkos.h index dac73aa7f68..bca63c65e25 100644 --- a/src/KOKKOS/pair_tip4p_kokkos.h +++ b/src/KOKKOS/pair_tip4p_kokkos.h @@ -389,7 +389,9 @@ class PairTIP4PKokkos : public PairCPUBase { { this->eflag = eflag_in; this->vflag = vflag_in; - this->ev_init(this->eflag,this->vflag); + // alloc = 0: the per-atom energy/virial arrays are allocated below through + // the Kokkos dual views, so Pair::ev_setup() must not allocate plain arrays + this->ev_init(this->eflag,this->vflag,0); this->atomKK->sync(this->execution_space,this->datamask_read); diff --git a/src/KOKKOS/region_block_kokkos.cpp b/src/KOKKOS/region_block_kokkos.cpp index 3802c22f968..4722797ed4e 100644 --- a/src/KOKKOS/region_block_kokkos.cpp +++ b/src/KOKKOS/region_block_kokkos.cpp @@ -16,7 +16,6 @@ #include "atom_kokkos.h" #include "atom_masks.h" -#include "memory_kokkos.h" using namespace LAMMPS_NS; using namespace MathSpecialKokkos; @@ -28,16 +27,6 @@ RegBlockKokkos::RegBlockKokkos(LAMMPS *lmp, int narg, char **arg) : RegBlock(lmp, narg, arg) { atomKK = (AtomKokkos*) atom; - memoryKK->create_kokkos(d_contact,6,"region_block:d_contact"); -} - -/* ---------------------------------------------------------------------- */ - -template -RegBlockKokkos::~RegBlockKokkos() -{ - if (copymode) return; - memoryKK->destroy_kokkos(d_contact); } /* ---------------------------------------------------------------------- */ @@ -47,6 +36,7 @@ void RegBlockKokkos::match_all_kokkos(int groupbit_in, DAT::tdual_in { groupbit = groupbit_in; d_match = k_match_in.template view(); + k_remap.setup(domain); auto execution_space = ExecutionSpaceFromDevice::space; atomKK->sync(execution_space, X_MASK | MASK_MASK); d_x = atomKK->k_x.view(); diff --git a/src/KOKKOS/region_block_kokkos.h b/src/KOKKOS/region_block_kokkos.h index 24285e2448f..c7767b574f8 100644 --- a/src/KOKKOS/region_block_kokkos.h +++ b/src/KOKKOS/region_block_kokkos.h @@ -28,6 +28,7 @@ RegionStyle(block/kk/host,RegBlockKokkos); #include "kokkos_base.h" #include "kokkos_type.h" #include "math_special_kokkos.h" +#include "region_remap_kokkos.h" namespace LAMMPS_NS { @@ -43,11 +44,21 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { typedef DeviceType device_type; typedef ArrayTypes AT; + // maximum number of contacts a particle can have with this region, i.e. the + // length of the contact list the caller has to provide. matches RegBlock::cmax + + static constexpr int MAXCONTACT = 6; + RegBlockKokkos(class LAMMPS *, int, char **); - ~RegBlockKokkos() override; void match_all_kokkos(int, DAT::tdual_int_1d) override; + void prematch() override + { + RegBlock::prematch(); + k_remap.setup(domain); + } + // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION void operator()(TagRegBlockMatchAll, const int&) const; @@ -56,6 +67,9 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { KOKKOS_INLINE_FUNCTION int match_kokkos(double x, double y, double z) const { + // Region::match() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x,y,z); if (dynamic) inverse_transform(x,y,z); if (openflag) return 1; return !(k_inside(x,y,z) ^ interior); @@ -63,12 +77,16 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE 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, Contact *contact) const { int ncontact; double xs, ys, zs; double xnear[3], xorig[3]; + // Region::surface() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x, y, z); + if (dynamic) { xorig[0] = x; xorig[1] = y; xorig[2] = z; inverse_transform(x, y, z); @@ -78,31 +96,34 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { if (!openflag) { if (interior) - ncontact = surface_interior_kokkos(xnear, cutoff); + ncontact = surface_interior_kokkos(xnear, cutoff, contact); else - ncontact = surface_exterior_kokkos(xnear, cutoff); + ncontact = surface_exterior_kokkos(xnear, cutoff, contact); } 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); + // most of the time, one of surface_int/ext() will return 0 + // however, when exactly on top of a periodic boundary + // both could return 1, so run exterior then interior + ncontact = surface_exterior_kokkos(xnear, cutoff, contact); + if (ncontact == 0) + ncontact = surface_interior_kokkos(xnear, cutoff, contact); } 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; + xs = xnear[0] - contact[i].delx; + ys = xnear[1] - contact[i].dely; + zs = xnear[2] - 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; + contact[i].delx = xorig[0] - xs; + contact[i].dely = xorig[1] - ys; + contact[i].delz = xorig[2] - zs; } } return ncontact; } - Kokkos::View d_contact; + RegionRemapKokkos k_remap; private: int groupbit; @@ -112,7 +133,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_interior_kokkos(double *x, double cutoff) + int surface_interior_kokkos(double *x, double cutoff, Contact *contact) const { double delta; @@ -126,58 +147,58 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { delta = x[0] - xlo; if (delta < cutoff && !open_faces[0]) { - d_contact[n].r = delta; - d_contact[n].delx = delta; - d_contact[n].dely = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 0; + contact[n].r = delta; + contact[n].delx = delta; + contact[n].dely = contact[n].delz = 0.0; + contact[n].radius = 0; + contact[n].iwall = 0; n++; } delta = xhi - x[0]; if (delta < cutoff && !open_faces[1]) { - d_contact[n].r = delta; - d_contact[n].delx = -delta; - d_contact[n].dely = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 1; + contact[n].r = delta; + contact[n].delx = -delta; + contact[n].dely = contact[n].delz = 0.0; + contact[n].radius = 0; + contact[n].iwall = 1; 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; + contact[n].r = delta; + contact[n].dely = delta; + contact[n].delx = contact[n].delz = 0.0; + contact[n].radius = 0; + 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; + contact[n].r = delta; + contact[n].dely = -delta; + contact[n].delx = contact[n].delz = 0.0; + contact[n].radius = 0; + 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; + contact[n].r = delta; + contact[n].delz = delta; + contact[n].delx = contact[n].dely = 0.0; + contact[n].radius = 0; + 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; + contact[n].r = delta; + contact[n].delz = -delta; + contact[n].delx = contact[n].dely = 0.0; + contact[n].radius = 0; + contact[n].iwall = 5; n++; } @@ -186,7 +207,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_exterior_kokkos(double *x, double cutoff) + int surface_exterior_kokkos(double *x, double cutoff, Contact *contact) const { double xp, yp, zp; double xc, yc, zc, dist, mindist; @@ -241,24 +262,24 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { if (mindist == MAXDOUBLEINT) return 0; } - add_contact(0, x, xp, yp, zp); - d_contact[0].iwall = 0; - if (d_contact[0].r < cutoff) return 1; + add_contact(0, x, xp, yp, zp, contact); + contact[0].iwall = 0; + if (contact[0].r < cutoff) return 1; return 0; } // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - void add_contact(int n, double *x, double xp, double yp, double zp) + void add_contact(int n, double *x, double xp, double yp, double zp, Contact *contact) const { double delx = x[0] - xp; double dely = x[1] - yp; 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; + contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); + contact[n].radius = 0; + contact[n].delx = delx; + contact[n].dely = dely; + contact[n].delz = delz; } // NOLINTNEXTLINE @@ -325,7 +346,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - void point_on_line_segment(double *a, double *b, double *c, double *d) + void point_on_line_segment(const double *a, const double *b, const double *c, double *d) const { double ba[3], ca[3]; @@ -349,7 +370,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - double inside_face(double *xproj, int iface) + double inside_face(double *xproj, int iface) const { if (iface < 2) { if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1; @@ -365,7 +386,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) + double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) const { double dot, d2, d2min; double xr[3], xproj[3], p[3]; diff --git a/src/KOKKOS/region_remap_kokkos.h b/src/KOKKOS/region_remap_kokkos.h new file mode 100644 index 00000000000..113295f5141 --- /dev/null +++ b/src/KOKKOS/region_remap_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. +------------------------------------------------------------------------- */ + +#ifndef LMP_REGION_REMAP_KOKKOS_H +#define LMP_REGION_REMAP_KOKKOS_H + +#include "domain.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +/* ---------------------------------------------------------------------- + device-side copy of Domain::remap(double *) + + Region::match() and Region::surface() map a coordinate back into the + periodic box before testing it, because not every region subclass handles + a region that reaches past a periodic edge. The KOKKOS regions run the + same test inside a kernel, where the Domain instance is not reachable, so + the few box quantities remap() needs are copied into this struct on the + host and captured by value with the region object. +------------------------------------------------------------------------- */ + +struct RegionRemapKokkos { + int triclinic, xperiodic, yperiodic, zperiodic; + double lo[3], hi[3], period[3]; + double boxlo[3], h[6], h_inv[6]; + + RegionRemapKokkos() : triclinic(0), xperiodic(0), yperiodic(0), zperiodic(0) + { + for (int i = 0; i < 3; i++) lo[i] = hi[i] = period[i] = boxlo[i] = 0.0; + for (int i = 0; i < 6; i++) h[i] = h_inv[i] = 0.0; + } + + // refresh from the Domain instance; call whenever the box may have changed + + void setup(Domain *domain) + { + triclinic = domain->triclinic; + xperiodic = domain->xperiodic; + yperiodic = domain->yperiodic; + zperiodic = domain->zperiodic; + for (int i = 0; i < 3; i++) { + if (triclinic) { + lo[i] = domain->boxlo_lamda[i]; + hi[i] = domain->boxhi_lamda[i]; + period[i] = domain->prd_lamda[i]; + } else { + lo[i] = domain->boxlo[i]; + hi[i] = domain->boxhi[i]; + period[i] = domain->prd[i]; + } + boxlo[i] = domain->boxlo[i]; + } + for (int i = 0; i < 6; i++) { + h[i] = domain->h[i]; + h_inv[i] = domain->h_inv[i]; + } + } + +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void remap(double &x, double &y, double &z) const + { + double coord[3]; + + if (triclinic == 0) { + coord[0] = x; coord[1] = y; coord[2] = z; + } else { + const double d0 = x - boxlo[0]; + const double d1 = y - boxlo[1]; + const double d2 = z - boxlo[2]; + coord[0] = h_inv[0]*d0 + h_inv[5]*d1 + h_inv[4]*d2; + coord[1] = h_inv[1]*d1 + h_inv[3]*d2; + coord[2] = h_inv[2]*d2; + } + + if (xperiodic) { + while (coord[0] < lo[0]) coord[0] += period[0]; + while (coord[0] >= hi[0]) coord[0] -= period[0]; + if (coord[0] < lo[0]) coord[0] = lo[0]; + } + + if (yperiodic) { + while (coord[1] < lo[1]) coord[1] += period[1]; + while (coord[1] >= hi[1]) coord[1] -= period[1]; + if (coord[1] < lo[1]) coord[1] = lo[1]; + } + + if (zperiodic) { + while (coord[2] < lo[2]) coord[2] += period[2]; + while (coord[2] >= hi[2]) coord[2] -= period[2]; + if (coord[2] < lo[2]) coord[2] = lo[2]; + } + + if (triclinic == 0) { + x = coord[0]; y = coord[1]; z = coord[2]; + } else { + x = h[0]*coord[0] + h[5]*coord[1] + h[4]*coord[2] + boxlo[0]; + y = h[1]*coord[1] + h[3]*coord[2] + boxlo[1]; + z = h[2]*coord[2] + boxlo[2]; + } + } +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/KOKKOS/region_sphere_kokkos.cpp b/src/KOKKOS/region_sphere_kokkos.cpp index b279891e1ba..a6b715dfaed 100644 --- a/src/KOKKOS/region_sphere_kokkos.cpp +++ b/src/KOKKOS/region_sphere_kokkos.cpp @@ -20,7 +20,6 @@ #include "atom_kokkos.h" #include "atom_masks.h" -#include "memory_kokkos.h" using namespace LAMMPS_NS; @@ -31,16 +30,6 @@ RegSphereKokkos::RegSphereKokkos(LAMMPS *lmp, int narg, char **arg) : RegSphere(lmp, narg, arg) { atomKK = (AtomKokkos*) atom; - memoryKK->create_kokkos(d_contact,1,"region_sphere:d_contact"); -} - -/* ---------------------------------------------------------------------- */ - -template -RegSphereKokkos::~RegSphereKokkos() -{ - if (copymode) return; - memoryKK->destroy_kokkos(d_contact); } /* ---------------------------------------------------------------------- */ @@ -50,6 +39,7 @@ void RegSphereKokkos::match_all_kokkos(int groupbit_in, DAT::tdual_i { groupbit = groupbit_in; d_match = k_match_in.template view(); + k_remap.setup(domain); auto execution_space = ExecutionSpaceFromDevice::space; atomKK->sync(execution_space, X_MASK | MASK_MASK); d_x = atomKK->k_x.view(); diff --git a/src/KOKKOS/region_sphere_kokkos.h b/src/KOKKOS/region_sphere_kokkos.h index bde0a9f6193..d36d694caca 100644 --- a/src/KOKKOS/region_sphere_kokkos.h +++ b/src/KOKKOS/region_sphere_kokkos.h @@ -27,6 +27,7 @@ RegionStyle(sphere/kk/host,RegSphereKokkos); #include "kokkos_base.h" #include "kokkos_type.h" +#include "region_remap_kokkos.h" namespace LAMMPS_NS { @@ -40,11 +41,21 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { typedef DeviceType device_type; typedef ArrayTypes AT; + // maximum number of contacts a particle can have with this region, i.e. the + // length of the contact list the caller has to provide. matches RegSphere::cmax + + static constexpr int MAXCONTACT = 1; + RegSphereKokkos(class LAMMPS *, int, char **); - ~RegSphereKokkos() override; void match_all_kokkos(int, DAT::tdual_int_1d) override; + void prematch() override + { + RegSphere::prematch(); + k_remap.setup(domain); + } + // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION void operator()(TagRegSphereMatchAll, const int&) const; @@ -53,6 +64,9 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { KOKKOS_INLINE_FUNCTION int match_kokkos(double x, double y, double z) const { + // Region::match() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x,y,z); if (dynamic) inverse_transform(x,y,z); if (openflag) return 1; return !(k_inside(x,y,z) ^ interior); @@ -60,44 +74,52 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE 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, Contact *contact) const { int ncontact; double xs, ys, zs; double xnear[3], xorig[3]; - xorig[0] = x; xorig[1] = y; xorig[2] = z; - if (dynamic) + // Region::surface() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x, y, z); + + if (dynamic) { + xorig[0] = x; xorig[1] = y; xorig[2] = 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); + if (interior) ncontact = surface_interior_kokkos(xnear, cutoff, contact); else - ncontact = surface_exterior_kokkos(xnear, cutoff); + ncontact = surface_exterior_kokkos(xnear, cutoff, contact); } 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); + // most of the time, one of surface_int/ext() will return 0 + // however, when exactly on top of a periodic boundary + // both could return 1, so run exterior then interior + ncontact = surface_exterior_kokkos(xnear, cutoff, contact); + if (ncontact == 0) + ncontact = surface_interior_kokkos(xnear, cutoff, contact); } 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; + xs = xnear[0] - contact[i].delx; + ys = xnear[1] - contact[i].dely; + zs = xnear[2] - 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; + contact[i].delx = xorig[0] - xs; + contact[i].dely = xorig[1] - ys; + contact[i].delz = xorig[2] - zs; } } return ncontact; } - Kokkos::View d_contact; + RegionRemapKokkos k_remap; private: int groupbit; @@ -107,7 +129,7 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_interior_kokkos(double *x, double cutoff) + int surface_interior_kokkos(double *x, double cutoff, Contact *contact) const { double delx = x[0] - xc; double dely = x[1] - yc; @@ -117,13 +139,13 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { double delta = radius - r; if (delta < cutoff) { - d_contact[0].r = delta; - d_contact[0].delx = delx * (1.0 - radius / r); - d_contact[0].dely = dely * (1.0 - radius / r); - d_contact[0].delz = delz * (1.0 - radius / r); - d_contact[0].radius = -radius; - d_contact[0].iwall = 0; - d_contact[0].varflag = 1; + contact[0].r = delta; + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); + contact[0].radius = -radius; + contact[0].iwall = 0; + contact[0].varflag = 1; return 1; } return 0; @@ -131,7 +153,7 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_exterior_kokkos(double *x, double cutoff) + int surface_exterior_kokkos(double *x, double cutoff, Contact *contact) const { double delx = x[0] - xc; double dely = x[1] - yc; @@ -141,13 +163,13 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { double delta = r - radius; if (delta < cutoff) { - d_contact[0].r = delta; - d_contact[0].delx = delx * (1.0 - radius / r); - d_contact[0].dely = dely * (1.0 - radius / r); - d_contact[0].delz = delz * (1.0 - radius / r); - d_contact[0].radius = radius; - d_contact[0].iwall = 0; - d_contact[0].varflag = 1; + contact[0].r = delta; + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); + contact[0].radius = radius; + contact[0].iwall = 0; + contact[0].varflag = 1; return 1; } return 0; @@ -155,16 +177,16 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - void add_contact(int n, double *x, double xp, double yp, double zp) + void add_contact(int n, double *x, double xp, double yp, double zp, Contact *contact) const { double delx = x[0] - xp; double dely = x[1] - yp; 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; + contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); + contact[n].radius = 0; + contact[n].delx = delx; + contact[n].dely = dely; + contact[n].delz = delz; } // NOLINTNEXTLINE diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index fbabd6c2bb0..736f866d664 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -62,6 +62,32 @@ struct Zero { } }; +/* ---------------------------------------------------------------------- + zero count entries of a per-atom array on both host sides, from first + + the loops that zero the device copy leave the host side alone. A style + that runs on the host adds its forces into the Kokkos host view, and a + style without KOKKOS support adds into the plain LAMMPS array behind it, + so whatever either of them added last step is still there and is counted + again when the two sides are brought together. Each call takes the same + range as the device loop it goes with, so that the atoms the plain code + leaves alone -- those outside an include group -- are left alone here too. +------------------------------------------------------------------------- */ + +template +static void zero_host_view(const View &v, int first, int count) +{ + Kokkos::parallel_for(Kokkos::RangePolicy(first,first+count),Zero(v)); +} + +template +static void zero_host(const DualView &k, int first, int count) +{ + if (count <= 0) return; + zero_host_view(k.view_hostkk(),first,count); + zero_host_view(k.view_host(),first,count); +} + /* ---------------------------------------------------------------------- */ VerletKokkos::VerletKokkos(LAMMPS *lmp, int narg, char **arg) : @@ -442,11 +468,13 @@ void VerletKokkos::run(int n) if (pair_compute_flag) { int prev_auto_sync = lmp->kokkos->auto_sync; if (!force->pair->kokkosable) lmp->kokkos->auto_sync = 1; - atomKK->sync(force->pair->execution_space,force->pair->datamask_read); + // the masked form only: the mask is the plain one when nothing is + // excluded, and syncing or claiming the full one first would put the + // force array back in play exactly where the overlap path is trying to + // keep it out atomKK->sync(force->pair->execution_space,~(~force->pair->datamask_read|datamask_exclude)); force->pair->compute(eflag,vflag); lmp->kokkos->auto_sync = prev_auto_sync; - atomKK->modified(force->pair->execution_space,force->pair->datamask_modify); atomKK->modified(force->pair->execution_space,~(~force->pair->datamask_modify|datamask_exclude)); timer->stamp(Timer::PAIR); } @@ -577,6 +605,16 @@ void VerletKokkos::force_clear() atomKK->k_f.clear_sync_state(); // ignore host forces/torques since device views atomKK->k_torque.clear_sync_state(); // will be cleared below + // the SPIN forces below are overwritten in the same way, so their host side + // has to be released here as well -- without this a host side left claimed + // from the setup is still claimed when the device side is claimed below, and + // the two disagree with nothing to say which one is current + + if (extraflag) { + atomKK->k_fm.clear_sync_state(); + atomKK->k_fm_long.clear_sync_state(); + } + // clear force on all particles // if either newton flag is set, also include ghosts // when using threads always clear all forces. @@ -586,10 +624,12 @@ void VerletKokkos::force_clear() if (force->newton) nall += atomKK->nghost; Kokkos::parallel_for(nall, Zero(atomKK->k_f.view_device())); + zero_host(atomKK->k_f,0,nall); atomKK->modified(Device,F_MASK); if (torqueflag) { Kokkos::parallel_for(nall, Zero(atomKK->k_torque.view_device())); + zero_host(atomKK->k_torque,0,nall); atomKK->modified(Device,TORQUE_MASK); } @@ -597,8 +637,10 @@ void VerletKokkos::force_clear() if (extraflag) { Kokkos::parallel_for(nall, Zero(atomKK->k_fm.view_device())); + zero_host(atomKK->k_fm,0,nall); atomKK->modified(Device,FM_MASK); Kokkos::parallel_for(nall, Zero(atomKK->k_fm_long.view_device())); + zero_host(atomKK->k_fm_long,0,nall); atomKK->modified(Device,FML_MASK); } @@ -608,10 +650,12 @@ void VerletKokkos::force_clear() } else { Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_f.view_device())); + zero_host(atomKK->k_f,0,atomKK->nfirst); atomKK->modified(Device,F_MASK); if (torqueflag) { Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_torque.view_device())); + zero_host(atomKK->k_torque,0,atomKK->nfirst); atomKK->modified(Device,TORQUE_MASK); } @@ -619,18 +663,22 @@ void VerletKokkos::force_clear() if (extraflag) { Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_fm.view_device())); + zero_host(atomKK->k_fm,0,atomKK->nfirst); atomKK->modified(Device,FM_MASK); Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_fm_long.view_device())); + zero_host(atomKK->k_fm_long,0,atomKK->nfirst); atomKK->modified(Device,FML_MASK); } if (force->newton) { auto range = Kokkos::RangePolicy(atomKK->nlocal, atomKK->nlocal + atomKK->nghost); Kokkos::parallel_for(range, Zero(atomKK->k_f.view_device())); + zero_host(atomKK->k_f,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,F_MASK); if (torqueflag) { Kokkos::parallel_for(range, Zero(atomKK->k_torque.view_device())); + zero_host(atomKK->k_torque,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,TORQUE_MASK); } @@ -638,8 +686,10 @@ void VerletKokkos::force_clear() if (extraflag) { Kokkos::parallel_for(range, Zero(atomKK->k_fm.view_device())); + zero_host(atomKK->k_fm,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,FM_MASK); Kokkos::parallel_for(range, Zero(atomKK->k_fm_long.view_device())); + zero_host(atomKK->k_fm_long,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,FML_MASK); } } diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index caacd4a297e..17b4fc5e2dd 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -155,6 +155,15 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : // register with Atom class reaxff = dynamic_cast(force->pair_match("^reaxff",0)); + + // this fix borrows the pair style's neighbor list when it can, to avoid + // building a second one. The KOKKOS version of pair reaxff builds its list + // on the device, where the plain ilist, numneigh and firstneigh below are + // not filled in, so fall back to the list this fix requests for itself -- + // the same one it uses when there is no reaxff pair style at all. + + if (reaxff && reaxff->kokkosable) reaxff = nullptr; + reaxflag = 0; nlevels_respa = 1; diff --git a/src/STUBS/mpi.h b/src/STUBS/mpi.h index 03b8d1fc263..d9a92dc778f 100644 --- a/src/STUBS/mpi.h +++ b/src/STUBS/mpi.h @@ -58,6 +58,7 @@ #define MPI_ANY_SOURCE (-1) #define MPI_STATUS_IGNORE NULL +#define MPI_STATUSES_IGNORE NULL #define MPI_Comm int #define MPI_Request int diff --git a/src/atom.h b/src/atom.h index b3bea4b1a6f..66633e2a2ff 100644 --- a/src/atom.h +++ b/src/atom.h @@ -356,10 +356,20 @@ class Atom : protected Pointers { void data_fix_compute_variable(int, int); virtual void allocate_type_arrays(); - void set_mass(const char *, int, const char *, int, int, int *); - void set_mass(const char *, int, int, double); - void set_mass(const char *, int, int, char **); - void set_mass(double *); + + // A style that reads or writes the plain per-atom arrays on the host says so + // with these, so that the KOKKOS package can bring that side up to date + // first and can carry the write over to the device afterwards. They do + // nothing without the package. Use them where a style touches the arrays + // directly and then calls something -- a force computation, the + // communication -- that works from the KOKKOS copies. + + virtual void sync_host_arrays(uint64_t) {} + virtual void modified_host_arrays(uint64_t) {} + virtual void set_mass(const char *, int, const char *, int, int, int *); + virtual void set_mass(const char *, int, int, double); + virtual void set_mass(const char *, int, int, char **); + virtual void set_mass(double *); void check_mass(const char *, int); int radius_consistency(int, double &); diff --git a/src/domain.h b/src/domain.h index 4b4835e2c2c..9cc54561f87 100644 --- a/src/domain.h +++ b/src/domain.h @@ -145,13 +145,13 @@ class Domain : protected Pointers { void closest_image(const double *const, const double *const, double *const); void remap(double *, imageint &); void remap(double *); - void remap_all(); + virtual void remap_all(); void remap_near(double *, double *); void unmap_inv(double *x, imageint); void unmap(double *, imageint); void unmap(const double *, imageint, double *); void unmap(const double *, const double *, imageint, int, double *, double *); - void image_flip(int, int, int); + virtual void image_flip(int, int, int); int ownatom(int, double *, imageint *, int); void define_general_triclinic(double *, double *, double *, double *); diff --git a/src/fix.cpp b/src/fix.cpp index 884ffa854be..e3cc98a6cca 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -203,7 +203,7 @@ void Fix::set_molecule(int, tagint, int, double *, double *, double *) energy is *much* faster. ------------------------------------------------------------------------- */ -void Fix::ev_setup(int eflag, int vflag) +void Fix::ev_setup(int eflag, int vflag, int alloc) { int i,n; @@ -234,18 +234,24 @@ void Fix::ev_setup(int eflag, int vflag) if (eflag_atom && atom->nlocal > maxeatom) { maxeatom = atom->nmax; - memory->destroy(eatom); - memory->create(eatom,maxeatom,"fix:eatom"); + if (alloc) { + memory->destroy(eatom); + memory->create(eatom,maxeatom,"fix:eatom"); + } } if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; - memory->destroy(vatom); - memory->create(vatom,maxvatom,6,"fix:vatom"); + if (alloc) { + memory->destroy(vatom); + memory->create(vatom,maxvatom,6,"fix:vatom"); + } } if (cvflag_atom && atom->nlocal > maxcvatom) { maxcvatom = atom->nmax; - memory->destroy(cvatom); - memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + if (alloc) { + memory->destroy(cvatom); + memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + } } // zero accumulators @@ -253,11 +259,11 @@ void Fix::ev_setup(int eflag, int vflag) // fixes tally it individually via fix_modify energy yes and compute_scalar() if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (eflag_atom) { + if (eflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) eatom[i] = 0.0; } - if (vflag_atom) { + if (vflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; @@ -268,7 +274,7 @@ void Fix::ev_setup(int eflag, int vflag) vatom[i][5] = 0.0; } } - if (cvflag_atom) { + if (cvflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { cvatom[i][0] = 0.0; @@ -291,7 +297,7 @@ void Fix::ev_setup(int eflag, int vflag) if thermo_virial is not set, virial tallying is disabled ------------------------------------------------------------------------- */ -void Fix::v_setup(int vflag) +void Fix::v_setup(int vflag, int alloc) { int i,n; @@ -310,19 +316,23 @@ void Fix::v_setup(int vflag) if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; - memory->destroy(vatom); - memory->create(vatom,maxvatom,6,"fix:vatom"); + if (alloc) { + memory->destroy(vatom); + memory->create(vatom,maxvatom,6,"fix:vatom"); + } } if (cvflag_atom && atom->nlocal > maxcvatom) { maxcvatom = atom->nmax; - memory->destroy(cvatom); - memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + if (alloc) { + memory->destroy(cvatom); + memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + } } // zero accumulators if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag_atom) { + if (vflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; @@ -333,7 +343,7 @@ void Fix::v_setup(int vflag) vatom[i][5] = 0.0; } } - if (cvflag_atom) { + if (cvflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { cvatom[i][0] = 0.0; diff --git a/src/fix.h b/src/fix.h index 5f259a666c3..ab4d6897b96 100644 --- a/src/fix.h +++ b/src/fix.h @@ -282,25 +282,34 @@ class Fix : protected Pointers { int dynamic; // recount atoms for temperature computes - void ev_init(int eflag, int vflag) + // alloc = 0 tells ev_setup()/v_setup() to only update the flags and the + // maxeatom/maxvatom/maxcvatom sizes, but not to allocate or zero any of the + // plain per-atom arrays. Styles that manage eatom/vatom themselves (e.g. + // the KOKKOS variants, which store them in dual views) must pass alloc = 0, + // since otherwise the plain arrays allocated here are orphaned by the style. + // Note that this covers cvatom as well: a style that passes alloc = 0 has to + // provide the centroid virial itself or set centroidstressflag to + // CENTROID_NOTAVAIL, since cvatom is then left unallocated. + + void ev_init(int eflag, int vflag, int alloc = 1) { if ((eflag && thermo_energy) || (vflag && thermo_virial)) - ev_setup(eflag, vflag); + ev_setup(eflag, vflag, alloc); else evflag = eflag_either = eflag_global = eflag_atom = eflag_only = vflag_either = vflag_global = vflag_atom = cvflag_atom = 0; } - void ev_setup(int, int); + void ev_setup(int, int, int alloc = 1); void ev_tally(int, int *, double, double, double *); - void v_init(int vflag) + void v_init(int vflag, int alloc = 1) { if (vflag && thermo_virial) - v_setup(vflag); + v_setup(vflag, alloc); else evflag = vflag_either = vflag_global = vflag_atom = cvflag_atom = 0; } - void v_setup(int); + void v_setup(int, int alloc = 1); void v_tally(int, int *, double, double *); void v_tally(int, int *, double, double *, int, int, int[][2], double *, double[][3]); void v_tally(int, int *, double, double *, double[][3], double[][3], double[]); diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index c95c57b9a44..259d283d4e8 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -778,6 +778,13 @@ void FixDeform::init() image flags to new values, making eqs in doc of Domain:image_flip incorrect ------------------------------------------------------------------------- */ +void FixDeform::migrate_atoms() +{ + irregular->migrate_atoms(); +} + +/* ---------------------------------------------------------------------- */ + void FixDeform::pre_exchange() { if (flip == 0) return; @@ -821,7 +828,7 @@ void FixDeform::pre_exchange() domain->remap_all(); domain->x2lamda(atom->nlocal); - irregular->migrate_atoms(); + migrate_atoms(); domain->lamda2x(atom->nlocal); flip = 0; diff --git a/src/fix_deform.h b/src/fix_deform.h index f36f4d85424..d79972ca75c 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -37,6 +37,10 @@ class FixDeform : public Fix { int setmask() override; void init() override; void pre_exchange() override; + + // the atom migration is a separate step so that an accelerator version can + // move the per-atom data to the side this runs on and back + virtual void migrate_atoms(); void end_of_step() override; void write_restart(FILE *) override; void restart(char *buf) override; diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index 24a5724cd69..655bf46949b 100644 --- a/src/fix_efield.cpp +++ b/src/fix_efield.cpp @@ -305,7 +305,7 @@ void FixEfield::post_force(int vflag) // reallocate efield array if necessary - if ((varflag == ATOM) && (atom->nmax > maxatom)) { + if (((varflag == ATOM) || (estyle == ATOM) || (pstyle == ATOM)) && (atom->nmax > maxatom)) { maxatom = atom->nmax; memory->destroy(efield); memory->create(efield, maxatom, 4, "efield:efield"); diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index e104714db88..47b0e791435 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -465,7 +465,7 @@ void FixWall::post_force(int vflag) { // virial setup - v_init(vflag); + v_setup_peratom(vflag); // energy intialize. // eflag is used to track whether wall energies have been communicated. diff --git a/src/fix_wall.h b/src/fix_wall.h index f863ccf182b..dd2b3aa84a9 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -47,6 +47,14 @@ class FixWall : public Fix { virtual void precompute(int) = 0; virtual void wall_particle(int, int, double) = 0; + + // set up the per-atom virial storage at the beginning of post_force(). + // accelerator styles that keep the per-atom virial in their own arrays + // override this and pass alloc = 0 to v_init(), so that the plain + // base-class vatom array is not allocated behind their back. + + virtual void v_setup_peratom(int vflag) { v_init(vflag); } + static void update_image_plane(int, int, double, double **, class Domain *); protected: diff --git a/src/respa.cpp b/src/respa.cpp index 0b08799d9b6..1e920f5891b 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -17,6 +17,7 @@ #include "respa.h" +#include "accelerator_kokkos.h" #include "angle.h" #include "atom.h" #include "atom_vec.h" @@ -288,6 +289,15 @@ void Respa::init() { Integrate::init(); + // rRESPA keeps its own copy of the forces of each level and clears and sums + // them through the plain LAMMPS arrays, and it calls the force computations + // without the transfers between the host and the device that run_style + // verlet does for the KOKKOS package. On a GPU that silently gives the + // wrong forces, so refuse the combination rather than run it. + + if (lmp->kokkos) + error->all(FLERR, "Run style respa is not supported by the KOKKOS package"); + // warn if no fixes if (modify->nfix == 0 && comm->me == 0) diff --git a/src/set.cpp b/src/set.cpp index 5fae968fb8d..4ca43ef5cac 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -16,6 +16,7 @@ #include "arg_info.h" #include "atom.h" +#include "atom_masks.h" #include "atom_vec.h" #include "atom_vec_body.h" #include "atom_vec_ellipsoid.h" @@ -534,6 +535,12 @@ void Set::selection(int n) void Set::invoke_actions() { + // every action below reads and writes the per-atom arrays through the plain + // pointers, so bring the host side up to date first and hand the writes over + // afterwards; without the KOKKOS package these do nothing + + atom->sync_host_arrays(ALL_MASK); + // reallocate per-atom variable storage if needed if (varflag && atom->nlocal > maxvariable) { @@ -593,6 +600,8 @@ void Set::invoke_actions() action->count_select = count_select; action->count_action = count_action; } + + atom->modified_host_arrays(ALL_MASK); } /* ---------------------------------------------------------------------- */