WIP: particel storage using Kokkos views - #5393
Open
RudolfWeeber wants to merge 164 commits into
Open
Conversation
Approved brainstorming outcome: replace the Particle struct (AoS) with a cell-sorted flat ParticleStore of component-major Kokkos View columns (per-field dual residency, ScatterView force/torque accumulation), with state, parameters, and observables in separate containers. Migration is incremental: fields are evicted from the struct group-by-group behind proxy accessors, with always-green tests and a 5% cumulative budget on LJ/P3M benchmarks. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- run_benchmarks: catch CalledProcessError from subprocess.run, print a clear one-line error naming the failed configuration, and return exit code 3 instead of letting the exception propagate (which would have aliased the compare regression exit code 1). - compare: append label= component to the printed configuration name when the label field is non-empty, so two groups differing only by label print distinct names. - module docstring: add exit-code table documenting codes 0/1/2/3. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The load-average check self-tripped: the gate's own 4-rank MPI configurations drive the 1-minute load average to ~4, which takes minutes to decay, so 'run' aborted after its first multi-rank configuration. Replace the /proc/loadavg reading with a two-sample measurement of 'foreign' CPU usage: utime+stime ticks from /proc/<pid>/stat summed over processes whose owner uid differs from os.getuid(), expressed as a percentage of one core. This ignores our own benchmark load and only reacts to other users on the shared machine. CLI flag --max-load becomes --max-foreign-cpu (default 50.0). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…torage Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…cleStorage Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- check_cell_storage_mutations.sh: rewrite as an embedded python3 scanner that reads each src/core file as one string and matches the mutation pattern across newlines, catching clang-format-wrapped calls that the old line-based grep missed. Document the alias limitation (mutation via a local reference is invisible; a green run is a tripwire, not a proof) and the decomposition swap/teardown exception. - ParticleListOperations.hpp: name the decomposition swap/teardown exception in the namespace doc (old cell storage is destroyed wholesale by destructors; phase 2 rebuilds the full store via mark-dirty). - benchmark_gate.py: 'run' refuses (exit 4) if the output CSV already exists, so stale rows cannot pollute the min-of-means comparison; clarify the exit-code table (exit 1 also covers missing configurations). - ParticleStore.hpp: spell accessor return types as std::size_t. - CellStructure.cpp: add the direct ParticleListOperations.hpp include. - ParticleListOperations_test.cpp: add extract-from-single-element and resize-ghost-storage-shrink test cases. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Replace reference bindings (auto &force = p.force()) with local copy + write-back, eliminate all force_and_torque() uses outside ghosts.cpp, and make binary-arithmetic uses of force()/torque() proxy-safe via explicit Utils::Vector3d conversions. Zero behavior change: today's accessors still return Utils::Vector3d&. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Bug 1: serialize_and_reduce FORCE branch gated on both policy==UPDATE and
direction==LOAD, so a plain MOVE+LOAD (e.g. lb_tracers' force ghost update)
fell into the else branch and serialized FROM the particle on an input
archive, discarding the received force. Fix by branching on direction first,
then on policy for the accumulate vs. assign distinction.
Bug 2: calc_transmit_size early-returned for FORCE with an assert that
data_parts==GHOSTTRANS_FORCE (false for mixed updates); in Release the assert
is a no-op and all non-FORCE parts were silently dropped from the size.
Replace with compositional sizing: compute force_size separately, mask FORCE
out of data_parts, run the Particle{}-based sizer only for remaining parts,
and return the sum.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add get_particle_force() and get_particle_torque_lab() in particle_node, each using a collective all_reduce over the live particle data on the owning rank. Route the Python-facing 'f' and 'torque_lab' getters in ParticleHandle through these new fetch functions instead of the fetch-cache detached Particle copy. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The ParticleForce member leaves the Particle struct in the same commit that makes the columns authoritative (spec section 4, single ownership). Non-const accessors return a write-through VectorReference; const accessors return values. Python getters fetch from the owning rank. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Phase-2 T7 moved force and torque out of the Particle struct into ParticleStore columns. Those Kokkos columns are rank-local and are not carried by Particle::serialize, which is used for whole-particle migration across MPI ranks during a global resort (AtomDecomposition and RegularDecomposition). A particle that migrated therefore arrived detached from the store, and the store rebuild zero-filled its new row, losing the force. This was invisible until a decomposition switch (global resort that moves particles) was followed by a force-reusing integrator step (run(0), which does not recompute forces): the reused force read 0. This surfaced as hybrid_decomposition::test_against_nsquare failing deterministically at 4 ranks (the n_square reference forces were wrong, not the hybrid ones); serial and regular-decomposition identity were unaffected because neither migrates a particle and then reuses forces. Fix: ferry force/torque with the migrating particle through a small transitional carrier on Particle (m_detached_force/m_detached_torque), serialized alongside the other fields; ParticleStore::assign_row seeds a detached (migrated or new) particle's row from the carrier. The column remains the source of truth once attached; brand-new particles keep the zero default. This restores the exact pre-migration behavior. The carrier is transitional and is removed in phase 7 when inter-rank exchange moves to per-field column packing. Adds a ParticleStore unit-test regression guard that serialize-roundtrips a particle and asserts assign_row seeds the migrated force/torque. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Final-review fixes for phase 2 of the ParticleStore migration: - VectorReference: add value-copying copy-assignment (write-through, like std::vector<bool>::reference); the copy constructor still rebinds. Add unit tests for value-copy and chained assignment. - particle_management test helper: reattach() re-points the moved particle at the store via attach_to_store instead of assign_row, so a real move no longer re-seeds from the stale migration carrier. - ShapeBasedConstraint: co-own the Kokkos runtime (shared_ptr<KokkosHandle>) captured in the lazy attach path, mirroring CellStructure, so the constraint's ParticleStore Views can be released before Kokkos::finalize even if it outlives the last CellStructure. - ParticleStore::assign_row: assert the row index is within bounds; record the adjudicated deviation in the design spec (id-based cross-check deferred to phase 5 when the id column exists). - Comment fixes: ParticleStore rebuild seeds new rows from the migration carrier (not zero-init); Particle_test clarifies force/torque ARE serialized via carriers but not applied to an already-attached particle. - Remove dead force_range() from ParticlePropertyIterator (no consumers). - trajectory_identity benchmark: clarify that P3M's required accuracy argument is unused with tune=False (kept: it is a required SI key; removing it breaks construction and the identity hashes are unchanged). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The H5MD "write" method reads local particle forces directly, but forces now live in the ParticleStore columns (migration phase 2) which are rebuilt lazily. Every other force-reading script-interface entry point already calls ensure_particle_store_synchronized(); the H5MD writer did not, so a write issued while the store is dirty (before the first integrator step, or right after a resort/particle addition with no intervening force read) could read stale/detached rows. Add the sync as the first step of the write branch. Also link the espresso_hdf5 target against Kokkos/Cabana: since force/torque moved into the ParticleStore, h5md_core.cpp transitively includes Kokkos_Core.hpp and no longer builds without those include paths. Add a Python regression test (testsuite/python/h5md.py) that marks the store dirty via a resort after adding a particle and then writes without any integrator step, checking the write does not crash and the fresh particle's force is zero. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ionReference (phase 3) Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add explicitly defaulted copy constructors to BasicVectorReference<T> and QuaternionReference next to their write-through copy assignments. The comment clarifies the asymmetry: the copy constructor rebinds the proxy (default pointer/stride copy), while copy assignment writes values through to the underlying storage. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…air kernel (simd ws1)" This reverts commit da29efd.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This reverts commit 0ec0035.
Restructure init_forces_and_thermostat into contiguous component sweeps under the component-major ParticleStore columns: - Force/torque init is now three (x,y,z) contiguous column copies (force = ext_force, torque = ext_torque) over [0, n_local), which auto-vectorize into packed AVX2 ymm moves, replacing the per-row Particle-proxy writes. External-force-disabled builds zero-fill the columns. The engine swim term (needs the quaternion-derived director) stays a per-row scalar update applied only when a swimmer exists. - The Langevin friction+noise pass is split off into a force-inlined per-row helper (apply_langevin_row). The id-keyed Philox draw stays scalar by design; per-particle arithmetic order is preserved (force = ext_force; force += (friction + noise)), so trajectories are bitwise-identical. always_inline (not [[gnu::flatten]]) is used because flattening this path fuses pref*v into the accumulate under -march=native and would break bitwise identity. Same-flags identity (native): lj/p3m/langevin all unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…simd ws3)" This reverts commit cccccc0.
This reverts commit f991360.
The z-line gather (simd ws2) wins single-threaded (-7.5% on the p3m benchmark at 1 rank, -3.8% at 4 MPI ranks, 100k particles) but measures slower when the assignment loop is OpenMP-threaded (+4.4% at 4 threads, 160k particles); multi-threaded execution keeps the per-point interpolate form. Same dispatch idiom as kokkos_parallel_range_for. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…lation units Move init_forces_and_thermostat (with its [[gnu::flatten]] Langevin tree) and the external_force helper out of forces.cpp into forces_init.cpp, and introduce a non-template update_verlet_state wrapper in short_range_verlet.cpp that is the sole instantiation site of update_cabana_state<VerletCriterion<>>. forces.cpp, energy.cpp and pressure.cpp now call the wrapper, so the AoSoA-commit / Verlet-list-build giant is deduped 3->1 and forces.cpp's per-TU inline-growth budget is no longer shared with the init tree. Pure code motion: bitwise identity holds on both canonical and native builds (lj + p3m). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Wave 1 of the short-range hot-path optimization: - Branchless cuboid minimum-image fold. Encode periodicity into a masked inverse box length (0 for non-periodic directions) so the per-component fold reduces to `dx - rint(dx * inv_masked) * L` with no branch. `rint` maps to a single rounding instruction. Results match the previous round-based fold across the pair-loop input domain (separations below 1.5 box lengths); confirmed bitwise-identical on the canonical lj and p3m trajectories. - Flat per-type-pair squared-cutoff table in VerletCriterion. The per-candidate cutoff query in the Verlet-list build becomes a dense table load instead of walking the InteractionsNonBonded pointer table; inactive pairs store a negative sentinel so the distance comparison rejects them without a separate activity check. - Hoist cuboid box parameters by value into the Verlet-build kernels via CuboidMinimumImage, instead of chasing the BoxGeometry reference for the box lengths on every candidate pair. - Force-inline Utils::Vector operator+= / operator-= (previously outlined as .isra clones called from inside the pair kernel). Canonical identity bitwise-preserved (lj, p3m); unit tests pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The `#pragma omp simd` on the fixed-trip charge-scatter and force-gather z-lines makes clang emit `-Wpass-failed=transform-warning` when it declines to vectorize the loop at the sanitizer build's `-O1`, which `-Werror` promotes to a hard failure. Restrict the pragma to non-clang compilers. The scatter loop has independent stores, so vectorization never affects its result. For the gather reduction, the pragma stays enabled on gcc (where it authorizes the z-line tree accumulation that fixes the canonical single-threaded gather); on clang the loop keeps the sequential order it already used, since the pragma was failing to apply there anyway. Results differ only at the rounding level, as the gather already does across compilers. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Wave 2 of the short-range hot-path optimization. Add SpecializedForcesKernel<HasCoulomb>, which owns the Verlet-list loop and runs once per particle instead of being invoked per pair by Cabana::neighbor_parallel_for. This lets it obtain the ScatterView accessor once per particle rather than once per pair -- the per-pair access() is an omp_get_thread_num call that the profile showed dominating the LJ pair cost -- and hoist the per-particle position, type and charge and the cuboid box parameters out of the neighbor loop. `if constexpr (HasCoulomb)` compiles the real-space electrostatics path in or out. forces.cpp installs it through a new optional ShortRangeVerletPairLoop hook on cabana_short_range, dispatched by create_specialized_verlet_pair_loop, which selects it only for a cuboid box with no NPT virial, dipolar or ELC kernel, DPD thermostat, Thole or Gay-Berne pair, or particle exclusion, and falls back to ForcesKernel otherwise. Two instantiations cover the hot cases: pure central radial forces and central-plus-real-space-electrostatics. The kernel keeps the generic per-pair scatter order (i before j, both scaled from the same pair force; no register accumulation), so results are bitwise-identical on a single thread. Canonical identity preserved (lj, p3m); unit tests pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Doxygen cannot resolve a \ref to the namespaced free function template detail::get_mi_coord_masked; render it as a code span instead. Comment-only. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Process each particle's Verlet neighbors in fixed-size tiles inside SpecializedForcesKernel: a scalar gather of the tile's neighbor positions into SoA scratch, a vectorized minimum-image pass (CuboidMinimumImage:: batch_vector_dist2) that folds the whole tile and squares the distances at once, then a scalar short-range-force pass that runs only for the pairs whose squared distance passes the cutoff gate. The vectorized pass compiles to packed AVX (the fold uses the branchless masked-inverse form). It yields the same per-pair fold vector and squared distance as the scalar path -- the squared distance accumulates from zero in component order, matching Utils::Vector::norm2 -- and the force pass keeps the same neighbor order and i-before-j scatter, so results are bitwise-identical. Re-evaluation of the earlier batched-minimum-image idea (which regressed as a whole-Verlet-list transform): batching per particle into small hot scratch and vectorizing only the distance gate, with the scalar force confined to the pairs that survive it, turns it into a measured speedup on lj and p3m. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add a ShortRangeOnly template parameter to VerletCriterion that drops the electrostatics / dipolar / collision early-accept branches from operator() via `if constexpr`. update_verlet_state now builds the criterion from the cutoffs (deriving skin from the cell structure) and selects the ShortRangeOnly variant when none of those cutoffs is active -- exactly the condition under which the removed branches would never have accepted a pair, so the per-candidate build result is unchanged while the hot build loop sheds the dead comparisons. The charge/dipole-moment loads guarded by those branches were already skipped at runtime by short-circuit evaluation (the effective cutoff is negative when the feature is inactive), so this removes dead branches rather than live loads; the measured effect is a small lj build-loop speedup, p3m unchanged. Canonical identity bitwise-preserved (lj, p3m); unit tests pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Doxygen cannot resolve a \ref to the class call operator from the class documentation; describe it in prose. Comment-only. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add an optional --mesh to p3m.py that fixes the (cubic) P3M mesh instead of tuning it; the remaining parameters (cao, r_cut, alpha) are still tuned for the prescribed mesh and the mesh tuning limits are dropped. 0 (default) keeps the auto-tuned behavior. Also relaxes the minimum-steps-per-tick assert from 50 to 20 so higher particle counts remain runnable. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
heFFTe 2.4.1 runs its single-rank FFT reshapes serially; on a single MPI rank that caps P3M's OpenMP scaling. Apply cmake/heffte.patch during the heFFTe FetchContent step to add `#pragma omp parallel for` to the pack/unpack and transpose reshape kernels and the post-FFT scaling in heffte_pack3d.h. Each output element is written exactly once from a fixed source, so the reshape is a pure permutation and the result is bitwise-unchanged. Measured (single rank, 10k particles, 1->4 OpenMP threads): p3m mesh=96 25.9->22.5 ms (efficiency 34%->40%), mesh=64 8.62->7.14 ms (41%->49%); 1-thread unchanged; trajectory identity bitwise-preserved. Patch mechanism mirrors the existing cmake/cabana.patch. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
extract_block and pad_with_zeros_discard_imag copy a mesh-sized field lane by lane on a single thread every force calculation, which the profile showed as a large serial `memmove` slice of P3M's single-rank cost. Replace the running destination/source iterators with explicit per-lane offsets so the outer loop runs under `#pragma omp parallel for`: each iteration touches a disjoint set of lanes, so the copy is thread-safe and the result is bitwise-identical (a pure permutation). Combined with the threaded heFFTe reshapes, single-rank P3M at 10k particles / 4 OpenMP threads improves further -- mesh=96 21.1->19.9 ms, mesh=64 7.3->6.7 ms; 1-thread unchanged; trajectory identity bitwise-preserved. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
extract_block writes every element of its output, so value-initializing the vector first is a wasted full-mesh memset each force calculation. Allocate it with a default-init allocator instead. Bitwise-identical. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Records the OpenMP thread-scaling investigation and its root cause (serial reciprocal-space path), the single-rank micro-optimizations (heFFTe reshape + crop/pad copy threading), and the MPI-rank strong-scaling confirmation that ranks — not threads — are the P3M parallelism lever (heFFTe domain-decomposes the mesh and transposes across ranks). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Add a kokkos-fft-based P3M FFT backend for single-MPI-rank runs, selected behind the CMake option ESPRESSO_BUILD_WITH_KOKKOS_FFT (default OFF). Multi-rank runs keep heFFTe (kokkos-fft is a local, non-distributed transform). Introduce an abstract P3MFFTBackend interface (the 6-method surface P3M already uses); heFFTe is exposed through a thin P3MFFTHeffte adapter around the unchanged P3MFFT, and P3MFFTKokkos implements the same interface with an r2c transform on Kokkos host views, reproducing heFFTe's row-major layout and its unscaled (scale::none) convention. CoulombP3MState::fft becomes a shared_ptr<P3MFFTBackend>; init_cpu_kernels picks the kokkos-fft backend when running CPU on a single rank. Forces match the heFFTe path to floating-point round-off (max rel diff 7e-16; coulomb_cloud_wall, p3m_madelung, p3m_fft pass) and the single-rank P3M step is ~7-10%% faster at mesh 64/96. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Execute the forward/backward transforms in place on the caller's buffers via kokkos-fft's new-array execute path instead of staging every transform through owned scratch views. Plans are built from the exact buffers they run on and cached by pointer, which both satisfies FFTW's alignment requirement (kokkos-fft plans with FFTW_ESTIMATE) and lets the stable k-space / no-halo real buffers be reused copy-free across steps. Only the transient forward input (a fresh extract_block buffer each step) is still staged into an aligned scratch. The c2r backward destroys its input, which is safe: ks_E_fields[d] is recomputed every step and not read afterwards. Forces still match heFFTe to round-off (max rel diff 7e-16; coulomb_cloud_wall, p3m_madelung pass); the single-rank P3M step is now ~14-22%% faster than heFFTe at mesh 64/96 (up from ~7-10%% with the staging copies). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Give the FFT backend a forward_input_buffer() and split extract_block into an extract_block_into(out, ...) primitive plus the allocating wrapper, so kernel_ks_charge_density extracts the no-halo density directly into the backend's own buffer. The kokkos-fft backend hands back its aligned scratch and then transforms in place, removing the last staging copy; the heFFTe adapter hands back a persistent buffer, dropping the per-step allocation. heFFTe output is bitwise identical; kokkos-fft still matches it to round-off (7e-16; coulomb_cloud_wall, p3m_madelung pass). Removing the copy trims another 0-5%% off the single-rank kokkos-fft P3M step (most at high thread counts). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Make the whole package compile and run under Kokkos_ENABLE_CUDA=ON with the core executing on the host, so the CUDA Kokkos backend and the legacy GPU actors coexist and the core->device migration can proceed incrementally. Every core Kokkos kernel dispatch and view is pinned to the host execution/memory space: DefaultExecutionSpace -> DefaultHostExecutionSpace; bare-count parallel_for/reduce given an explicit host RangePolicy; the kokkos_parallel_range_for default policy and the for_each_3d MDRangePolicy pinned to host; and the force/virial/id-index/bond/energy/pressure Views (and ScatterViews) given explicit Kokkos::HostSpace -- they previously defaulted to DefaultMemorySpace = CudaSpace under CUDA and faulted when host code touched them. On the CPU build DefaultHostExecutionSpace/HostSpace resolve to the same OpenMP/HostSpace types, so this is a no-op there: lj and p3m trajectory identity hashes are unchanged (1671c333/6ac402d2, d2dfa428/e83aba7b at 1/4 threads). Under Kokkos_ENABLE_CUDA=ON (clang-19, sm_86) the package now builds clean and an LJ integration runs correctly on the host core. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
… phase 1) Add device-space view accessors for the hot state columns (position, velocity, force, image_box, id, type, q, mass) and explicit sync_state_to_device / sync_state_to_host helpers that deep_copy them between the DualView host and device mirrors. This is the data-residency substrate the device core kernels need; it is inert on the host path (the syncs run only on the opt-in device execution path, and on a host-only Kokkos build the device mirror aliases the host mirror). CPU identity unchanged (lj 1671c333, p3m d2dfa428). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ore phase 3) First Kokkos device kernel of the SoA core. A device short-range LJ pair force runs on the GPU when ESPRESSO_GPU_CORE=1 and the system is plain LJ on a cuboid box (no Coulomb/dipoles/ELC/DPD/Thole/Gay-Berne, no exclusions, no NPT virial, single rank); otherwise the existing host path is used unchanged. Implementation is isolated in its own CUDA translation unit (forces_lj_device.cu, compiled by clang-CUDA) and declared to forces.cpp via short_range_cabana.hpp, so forces.cpp gains no new include -- pulling the Kokkos/Cabana headers into that TU perturbed its Coulomb short-range FP codegen and broke P3M bitwise identity; keeping the device code in a separate TU keeps forces.cpp byte-identical. Hooks in through the existing ShortRangeVerletPairLoop std::function. A flat device LJ param table + a device-clean pair kernel (double precision, atomic half-list scatter, self-contained min-image fold, box params host-computed) write into local_force; bonds stay in scatter_force so the existing reduce sums them. Validated (clang-19, sm_86): device build clean; device vs host LJ forces match to 6e-14 (atomic-order only); runs on the GPU (nvidia-smi). CPU path unchanged and bitwise-identical (lj 1671c333, p3m d2dfa428). File/factory named for the full short-range kernel; LJ is the first implemented potential. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…re phase 3.1) Extend the opt-in device short-range pair force to add the P3M real-space Coulomb contribution, so a P3M simulation runs its short-range (real space) on the GPU alongside LJ. A device-callable coulomb factor reproduces CoulombP3M::pair_force (USE_ERFC_APPROXIMATION=1 branch, via the constexpr Utils::AS_erfc_part) using the already-synced device q column; the gate now accepts a P3M Coulomb solver (extracting prefactor/alpha/r_cut) and still falls back to host for any non-P3M Coulomb, ELC, dipoles, DPD, Thole/Gay-Berne, exclusions, NPT or non-cuboid box. LJ and Coulomb each self-gate on their own cutoff and share the neighbor list. Change is confined to forces_lj_device.cu (the isolated device TU), so forces.cpp and the CPU build are untouched (identity preserved). Validated (clang-19, sm_86): device vs host LJ+P3M-real-space forces match to 2e-14 (atomic-order only). Perf is not yet a win (double precision on a consumer GPU + per-step transfers); single precision and transfer reuse are later increments. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The opt-in device (Kokkos CUDA) short-range pair loop allocated five device Views and synced all eight ParticleStore columns to the device on every force call, leaving it ~2x slower than the host path even though the kernel itself is fast. Per-step allocation and over-copy were the bottleneck, not the arithmetic (single precision did not help). Keep the Verlet list build on the host and copy only what changes: - Hold persistent device buffers (counts/neighbors/type/row_map/force) in an opaque DeviceShortRangeBuffers owned via a shared_ptr on CellStructure, reused across steps and reallocated only when a dimension changes. - Copy the neighbor list and pack maps to the device only on a Verlet rebuild, keyed on a new monotonic Cabana Verlet-list generation counter (bumped in rebuild_verlet_list_cabana). - Sync only the per-step columns the kernel reads: positions every step, charges only when Coulomb is active (was: all eight columns). - Re-zero the persistent force accumulator in place each step. The buffers type is forward-declared in CellStructure.hpp and defined in the .cu, so no Kokkos headers reach the shared header (that perturbs host FP codegen). The shared_ptr's type-erased deleter keeps CellStructure.cpp free of the complete type; CellStructure's KokkosHandle guarantees the buffers are destroyed before Kokkos::finalize. Device vs host (RTX 3070 Ti, double precision, ms/step): LJ N=1728 : 0.22 vs 0.14 (small-N, kernel-launch bound) P3M N=1728 : 0.74 vs 0.83 (1.12x faster) LJ N=10648 : 0.85 vs 0.92 (1.08x faster) LJ N=21952 : 1.67 vs 1.87 (1.12x faster) Previously the device path was ~2x slower across the board. Device forces stay bit-for-bit close to the host: max abs 5e-16 single step, 3.5e-15 over 150 steps across many Verlet rebuilds (exercising buffer reuse and copy-on-rebuild). CPU path unchanged: LJ/P3M trajectory hashes identical (1671c333, d2dfa428). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This is a highly experimental agent run push to switch the particle storage to Kokkos views entirely, eliminating the Particle struct in the core.
Python interface stays as it is.
We use Fable for planning/design (while we have it) and Opus 4.8 for the actual coding tassks.