From c66f0e3fe98571053008d653e15d69e0d38e7c42 Mon Sep 17 00:00:00 2001 From: stekajack Date: Thu, 8 May 2025 12:05:26 +0200 Subject: [PATCH 01/16] WIP: n_replicas added to DDS_GPU --- .../magnetostatics/dipolar_direct_sum_gpu.cpp | 8 +- .../magnetostatics/dipolar_direct_sum_gpu.hpp | 3 +- .../dipolar_direct_sum_gpu_cuda.cu | 131 +++++++++++++----- .../dipolar_direct_sum_gpu_cuda.cuh | 8 +- src/python/espressomd/magnetostatics.py | 2 +- .../magnetostatics/DipolarDirectSumGpu.hpp | 10 +- 6 files changed, 122 insertions(+), 40 deletions(-) diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index 63f99e24ca3..2f695e8e89c 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -37,8 +37,12 @@ static void get_simulation_box(BoxGeometry const &box_geo, float *box, } } -DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor) { +DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor, int n_replicas) { set_prefactor(prefactor); + this->n_replicas = n_replicas; + if (n_replicas < 0) { + throw std::domain_error("Parameter 'n_replicas' must be >= 0"); + } } void DipolarDirectSumGpu::on_activation() const { @@ -66,7 +70,7 @@ void DipolarDirectSumGpu::add_long_range_forces() const { auto const dipoles_device = gpu.get_particle_dipoles_device(); DipolarDirectSum_kernel_wrapper_force( static_cast(prefactor), npart, positions_device, dipoles_device, - forces_device, torques_device, box, periodicity); + forces_device, torques_device, box, periodicity, n_replicas); } void DipolarDirectSumGpu::long_range_energy() const { diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.hpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.hpp index 2e05b702673..45622bd8ecb 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.hpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.hpp @@ -26,7 +26,8 @@ #include "magnetostatics/actor.hpp" struct DipolarDirectSumGpu : public Dipoles::Actor { - DipolarDirectSumGpu(double prefactor); + int n_replicas; + DipolarDirectSumGpu(double prefactor, int n_replicas); void on_activation() const; void on_boxl_change() const {} diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu index bff0f93c323..7574be5209e 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu @@ -56,14 +56,9 @@ __device__ inline void get_mi_vector_dds(float res[3], float const a[3], } } -__device__ void dipole_ia_force(float pf, float const *r1, float const *r2, - float const *dip1, float const *dip2, float *f1, - float *torque1, float *torque2, float box_l[3], - int periodic[3]) { - // Distance between particles - float dr[3]; - get_mi_vector_dds(dr, r1, r2, box_l, periodic); - +__device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1, + float const *dip2, float *f1, float *torque1, + float *torque2) { // Powers of distance auto const r_sq = scalar_product(dr, dr); auto const r_sq_inv = 1.0f / r_sq; @@ -134,20 +129,21 @@ __device__ float dipole_ia_energy(float pf, float const *r1, float const *r2, } // LCOV_EXCL_STOP +// constant‐memory shift list +__constant__ Int3 d_imageShifts[4096]; +__constant__ int d_nShifts; __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, float *pos, float *dip, float *f, - float *torque, float box_l[3], - int periodic[3]) { - + float *torque, + float const box_l[3], + int const periodic[3]) { auto const i = blockIdx.x * blockDim.x + threadIdx.x; - if (i >= n) return; // Kahan summation based on the wikipedia article // Force float fi[3], fsum[3], tj[3]; - // Torque float ti[3], tsum[3]; @@ -159,32 +155,74 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, // to global memory at the end. // Clear summation vars - for (unsigned int j = 0; j < 3; j++) { - // Force - fsum[j] = 0; - // Torque - tsum[j] = 0; + for (unsigned int k = 0; k < 3; ++k) { + fsum[k] = 0.0f; + tsum[k] = 0.0f; } - for (unsigned int j = i + 1; j < n; j++) { - dipole_ia_force(pf, pos + 3 * i, pos + 3 * j, dip + 3 * i, dip + 3 * j, fi, - ti, tj, box_l, periodic); - for (unsigned int k = 0; k < 3; k++) { - // Add rhs to global memory - atomicAdd(f + 3 * j + k, -fi[k]); - atomicAdd((torque + 3 * j + k), tj[k]); - tsum[k] += ti[k]; + // --- Self‐images of particle i (all shifts except the primary) --- + // these only update thread‐private fsum/tsum, so no atomics + for (int s = 1; s < d_nShifts; ++s) { + Int3 sh = d_imageShifts[s]; + float dr[3] = {sh.x * box_l[0], sh.y * box_l[1], sh.z * box_l[2]}; + dipole_ia_force(pf, dr, dip + 3 * i, dip + 3 * i, fi, ti, tj); + for (int k = 0; k < 3; ++k) { fsum[k] += fi[k]; + tsum[k] += ti[k]; } } - // Add the left hand side result to global memory - for (int j = 0; j < 3; j++) { - atomicAdd(f + 3 * i + j, fsum[j]); - atomicAdd(torque + 3 * i + j, tsum[j]); + // Pre‐load pos[i] once + float pi[3] = {pos[3 * i], pos[3 * i + 1], pos[3 * i + 2]}; + + // --- Pairwise interactions i ↔ j, for j>i, summing over all images --- + for (unsigned int j = i + 1; j < n; ++j) { + // Base MIC vector in the primary box + float pj[3] = {pos[3 * j], pos[3 * j + 1], pos[3 * j + 2]}; + float d0[3]; + get_mi_vector_dds(d0, pi, pj, box_l, periodic); + + // Loop over every image shift + for (int s = 0; s < d_nShifts; ++s) { + Int3 sh = d_imageShifts[s]; + float dr[3]; + if (sh.x == 0 && sh.y == 0 && sh.z == 0) { + // zero‐shift: use true minimal‐image vector + dr[0] = d0[0]; + dr[1] = d0[1]; + dr[2] = d0[2]; + } else { + // explicit replica: offset the MIC vector + dr[0] = d0[0] + sh.x * box_l[0]; + dr[1] = d0[1] + sh.y * box_l[1]; + dr[2] = d0[2] + sh.z * box_l[2]; + } + + // Compute interactions for this image + dipole_ia_force(pf, dr, dip + 3 * i, dip + 3 * j, fi, ti, tj); + + // 1) deposit force/torque on particle j **atomically** + // (many threads i racing to update the same j) + for (int k = 0; k < 3; ++k) { + atomicAdd(f + 3 * j + k, -fi[k]); + atomicAdd(torque + 3 * j + k, tj[k]); + } + + // 2) accumulate into thread‐private sums for particle i + for (int k = 0; k < 3; ++k) { + fsum[k] += fi[k]; + tsum[k] += ti[k]; + } + } } -} + // Add i’s total to global memory **atomically**, + // in case any other asynchronous update might race + for (int k = 0; k < 3; ++k) { + atomicAdd(f + 3 * i + k, fsum[k]); + atomicAdd(torque + 3 * i + k, tsum[k]); + } +} // LCOV_EXCL_START __device__ void dds_sumReduction(float *input, float *sum) { auto const tid = static_cast(threadIdx.x); @@ -244,17 +282,44 @@ inline void copy_box_data(float **box_l_gpu, int **periodic_gpu, cudaMemcpy(*periodic_gpu, periodic, s_per, cudaMemcpyHostToDevice)); } +// build all shifts with x²+y²+z² ≤ n_replicas² +static std::vector build_shifts(int n_replicas, int periodic[3]) { + std::vector shifts; + int rx = n_replicas * periodic[0]; + int ry = n_replicas * periodic[1]; + int rz = n_replicas * periodic[2]; + int cutoff2 = n_replicas * n_replicas; + shifts.reserve((2 * rx + 1) * (2 * ry + 1) * (2 * rz + 1)); + for (int x = -rx; x <= rx; ++x) + for (int y = -ry; y <= ry; ++y) + for (int z = -rz; z <= rz; ++z) + if (x * x + y * y + z * z <= cutoff2) + shifts.push_back({x, y, z}); + return shifts; +} + void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, float *dip, float *f, float *torque, - float box_l[3], int periodic[3]) { + float box_l[3], int periodic[3], + int n_replicas) { - unsigned int const bs = 64; + unsigned int const bs = 32; dim3 grid(1, 1, 1); dim3 block(1, 1, 1); if (n == 0) return; + static bool first = true; + if (first) { + auto shifts = build_shifts(n_replicas, periodic); + cudaMemcpyToSymbol(d_imageShifts, shifts.data(), + shifts.size() * sizeof(Int3)); + int nShifts = static_cast(shifts.size()); + cudaMemcpyToSymbol(d_nShifts, &nShifts, sizeof(int)); + first = false; + } + if (n <= bs) { grid.x = 1; block.x = n; diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh index f238f954eae..8545927328c 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh @@ -23,11 +23,17 @@ #ifdef DIPOLAR_DIRECT_SUM +// simple 3-int struct for image shifts +struct Int3 { + int x, y, z; +}; + void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, float *pos, float *dip, float box_l[3], int periodic[3], float *E); void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, float *dip, float *f, float *torque, - float box_l[3], int periodic[3]); + float box_l[3], int periodic[3], + int n_replicas); #endif // DIPOLAR_DIRECT_SUM diff --git a/src/python/espressomd/magnetostatics.py b/src/python/espressomd/magnetostatics.py index 1884ea884c7..8260a3f6fa1 100644 --- a/src/python/espressomd/magnetostatics.py +++ b/src/python/espressomd/magnetostatics.py @@ -255,7 +255,7 @@ class DipolarDirectSumGpu(MagnetostaticInteraction): _so_features = ("DIPOLAR_DIRECT_SUM", "CUDA") def default_params(self): - return {} + return {"n_replicas": 0} def required_keys(self): return {"prefactor"} diff --git a/src/script_interface/magnetostatics/DipolarDirectSumGpu.hpp b/src/script_interface/magnetostatics/DipolarDirectSumGpu.hpp index fca8bb74949..d197483d38a 100644 --- a/src/script_interface/magnetostatics/DipolarDirectSumGpu.hpp +++ b/src/script_interface/magnetostatics/DipolarDirectSumGpu.hpp @@ -38,12 +38,18 @@ namespace Dipoles { class DipolarDirectSumGpu : public Actor { public: - DipolarDirectSumGpu() = default; + DipolarDirectSumGpu() { + add_parameters({ + {"n_replicas", AutoParameter::read_only, + [this]() { return actor()->n_replicas; }}, + }); + } void do_construct(VariantMap const ¶ms) override { context()->parallel_try_catch([this, ¶ms]() { m_actor = std::make_shared( - get_value(params, "prefactor")); + get_value(params, "prefactor"), + get_value(params, "n_replicas")); }); } }; From f72a175fc0f3d24e30d1a2a1d4c9f13a0f775817 Mon Sep 17 00:00:00 2001 From: stekajack Date: Thu, 8 May 2025 15:17:43 +0000 Subject: [PATCH 02/16] test for rank consistency dds GPU --- testsuite/python/dipolar_direct_summation.py | 36 ++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/testsuite/python/dipolar_direct_summation.py b/testsuite/python/dipolar_direct_summation.py index 61b3af5f48a..cceec180717 100644 --- a/testsuite/python/dipolar_direct_summation.py +++ b/testsuite/python/dipolar_direct_summation.py @@ -296,6 +296,42 @@ def test_inner_loop_consistency_cpu(self): np.testing.assert_allclose(node_01_energy, node_00_energy, **tol) np.testing.assert_allclose(node_01_forces, node_00_forces, **tol) np.testing.assert_allclose(node_01_torques, node_00_torques, **tol) + + def test_inner_loop_consistency_gpu(self): + system = self.system + system.periodicity = [True, True, True] + tol = {"atol": 1e-6, "rtol": 1e-6} + p1 = system.part.add(pos=[0., 0., 0.], dip=[0., 0., 1.], + rotation=[True, True, True]) + p2 = system.part.add(pos=[1., 0., 0.], dip=[0., 0., 1.], + rotation=[True, True, True]) + for n_replicas in [0, 1]: + system.magnetostatics.clear() + solver = espressomd.magnetostatics.DipolarDirectSumGpu( + prefactor=1., n_replicas=n_replicas) + system.magnetostatics.solver = solver + + # intra-node calculation + p1.pos = [system.box_l[0] / 2. - 0.1, 0., 2.] + p2.pos = [system.box_l[0] / 2. + 0.1, 0., 0.] + system.integrator.run(steps=0, recalc_forces=True) + assert p1.node != p2.node + node_01_energy = system.analysis.energy()["dipolar"] + node_01_forces = np.copy(system.part.all().f) + node_01_torques = np.copy(system.part.all().torque_lab) + + # inter-node calculation + p1.pos = [0.1, 0., 2.] + p2.pos = [0.3, 0., 0.] + system.integrator.run(steps=0, recalc_forces=True) + assert p1.node == p2.node + node_00_energy = system.analysis.energy()["dipolar"] + node_00_forces = np.copy(system.part.all().f) + node_00_torques = np.copy(system.part.all().torque_lab) + + np.testing.assert_allclose(node_01_energy, node_00_energy, **tol) + np.testing.assert_allclose(node_01_forces, node_00_forces, **tol) + np.testing.assert_allclose(node_01_torques, node_00_torques, **tol) if __name__ == "__main__": From b605624628c81500a34a14cb326e9fc6fd4cb011 Mon Sep 17 00:00:00 2001 From: stekajack Date: Thu, 22 May 2025 11:48:18 +0000 Subject: [PATCH 03/16] dip_fld_trackiing added --- src/core/forces.cpp | 18 ++ .../magnetostatics/dipolar_direct_sum_gpu.cpp | 13 +- .../dipolar_direct_sum_gpu_cuda.cu | 257 +++++++++++------- .../dipolar_direct_sum_gpu_cuda.cuh | 8 +- src/core/magnetostatics/dipoles.cpp | 2 + src/core/magnetostatics/dipoles.hpp | 2 + src/core/system/GpuParticleData.cpp | 34 +++ src/core/system/GpuParticleData.hpp | 13 +- src/core/system/GpuParticleData_cuda.cu | 57 ++++ testsuite/python/dipole_field_tracking.py | 8 + 10 files changed, 309 insertions(+), 103 deletions(-) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 148eaf8cd84..a156ec42b86 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -97,6 +97,16 @@ void init_forces(const CellStructure &cell_structure) { init_forces_ghosts(cell_structure); } +#ifdef DIPOLE_FIELD_TRACKING +void invalidate_dip_fld(const CellStructure &cell_structure) { +#ifdef CALIPER + CALI_CXX_MARK_FUNCTION; +#endif + + cell_structure.for_each_local_particle( + [](Particle &p) { p.dip_fld() = {0., 0., 0.}; }); +} +#endif void init_forces_ghosts(const CellStructure &cell_structure) { cell_structure.for_each_ghost_particle( @@ -151,6 +161,11 @@ void System::System::calculate_forces() { } #endif init_forces(*cell_structure); +#ifdef DIPOLE_FIELD_TRACKING + // reset dipole field + invalidate_dip_fld(*cell_structure); +#endif + thermostat_force_init(); calc_long_range_forces(particles); @@ -224,6 +239,9 @@ void System::System::calculate_forces() { CALI_MARK_BEGIN("copy_forces_from_GPU"); #endif gpu.copy_forces_to_host(particles, this_node); +#ifdef DIPOLE_FIELD_TRACKING + gpu.copy_dip_fld_to_host(particles, this_node); +#endif #ifdef CALIPER CALI_MARK_END("copy_forces_from_GPU"); #endif diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index 2f695e8e89c..420a6e21fed 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -51,6 +51,9 @@ void DipolarDirectSumGpu::on_activation() const { gpu_particle_data.enable_property(GpuParticleData::prop::torque); gpu_particle_data.enable_property(GpuParticleData::prop::pos); gpu_particle_data.enable_property(GpuParticleData::prop::dip); +#ifdef DIPOLE_FIELD_TRACKING + gpu_particle_data.enable_property(GpuParticleData::prop::dip_fld); +#endif } void DipolarDirectSumGpu::add_long_range_forces() const { @@ -68,8 +71,16 @@ void DipolarDirectSumGpu::add_long_range_forces() const { auto const torques_device = gpu.get_particle_torques_device(); auto const positions_device = gpu.get_particle_positions_device(); auto const dipoles_device = gpu.get_particle_dipoles_device(); +#ifdef DIPOLE_FIELD_TRACKING + auto const dipole_fields_device = gpu.get_particle_dip_fld_device(); +#endif DipolarDirectSum_kernel_wrapper_force( - static_cast(prefactor), npart, positions_device, dipoles_device, + static_cast(prefactor), npart, positions_device, dipoles_device +#ifdef DIPOLE_FIELD_TRACKING + , + dipole_fields_device +#endif + , forces_device, torques_device, box, periodicity, n_replicas); } diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu index 7574be5209e..dab708f3059 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu @@ -58,7 +58,12 @@ __device__ inline void get_mi_vector_dds(float res[3], float const a[3], __device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1, float const *dip2, float *f1, float *torque1, - float *torque2) { + float *torque2 +#ifdef DIPOLE_FIELD_TRACKING + , + float *dip_fld1, float *dip_fld2 +#endif +) { // Powers of distance auto const r_sq = scalar_product(dr, dr); auto const r_sq_inv = 1.0f / r_sq; @@ -73,6 +78,18 @@ __device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1, auto const pe3 = scalar_product(dip2, dr); auto const pe4 = 3.0f * r5_inv; +#ifdef DIPOLE_FIELD_TRACKING + auto const rep1 = pe3 * pe4; + auto const rep2 = pe2 * pe4; + dip_fld1[0] = pf * (rep1 * dr[0] - dip2[0] * r3_inv); + dip_fld1[1] = pf * (rep1 * dr[1] - dip2[1] * r3_inv); + dip_fld1[2] = pf * (rep1 * dr[2] - dip2[2] * r3_inv); + + dip_fld2[0] = pf * (rep2 * dr[0] - dip1[0] * r3_inv); + dip_fld2[1] = pf * (rep2 * dr[1] - dip1[1] * r3_inv); + dip_fld2[2] = pf * (rep2 * dr[2] - dip1[2] * r3_inv); +#endif + // Force auto const aa = pe4 * pe1; auto const bb = -15.0f * pe2 * pe3 * r7_inv; @@ -129,89 +146,140 @@ __device__ float dipole_ia_energy(float pf, float const *r1, float const *r2, } // LCOV_EXCL_STOP -// constant‐memory shift list -__constant__ Int3 d_imageShifts[4096]; -__constant__ int d_nShifts; __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, - float *pos, float *dip, float *f, - float *torque, - float const box_l[3], - int const periodic[3]) { - auto const i = blockIdx.x * blockDim.x + threadIdx.x; + const float *pos, const float *dip +#ifdef DIPOLE_FIELD_TRACKING + , + float *dip_fld +#endif + , + float *f, float *torque, + const float box_l[3], + const int periodic[3], + int n_replicas) { + unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= n) return; - // Kahan summation based on the wikipedia article - // Force - float fi[3], fsum[3], tj[3]; - // Torque - float ti[3], tsum[3]; - - // There is one thread per particle. Each thread computes interactions - // with particles whose id is smaller than the thread id. - // The force and torque of all the interaction partners of the current thread - // is atomically added to global results at once. - // The result for the particle id equal to the thread id is atomically added - // to global memory at the end. - - // Clear summation vars - for (unsigned int k = 0; k < 3; ++k) { - fsum[k] = 0.0f; - tsum[k] = 0.0f; - } + // Load particle i + float xi = pos[3 * i + 0], yi = pos[3 * i + 1], zi = pos[3 * i + 2]; - // --- Self‐images of particle i (all shifts except the primary) --- - // these only update thread‐private fsum/tsum, so no atomics - for (int s = 1; s < d_nShifts; ++s) { - Int3 sh = d_imageShifts[s]; - float dr[3] = {sh.x * box_l[0], sh.y * box_l[1], sh.z * box_l[2]}; - dipole_ia_force(pf, dr, dip + 3 * i, dip + 3 * i, fi, ti, tj); - for (int k = 0; k < 3; ++k) { - fsum[k] += fi[k]; - tsum[k] += ti[k]; - } - } + const float *mi = dip + 3 * i; - // Pre‐load pos[i] once - float pi[3] = {pos[3 * i], pos[3 * i + 1], pos[3 * i + 2]}; + // Per‐thread accumulators + float fsum[3] = {0.0f, 0.0f, 0.0f}; + float tsum[3] = {0.0f, 0.0f, 0.0f}; +#ifdef DIPOLE_FIELD_TRACKING + float dfsum[3] = {0.0f, 0.0f, 0.0f}; +#endif - // --- Pairwise interactions i ↔ j, for j>i, summing over all images --- - for (unsigned int j = i + 1; j < n; ++j) { - // Base MIC vector in the primary box - float pj[3] = {pos[3 * j], pos[3 * j + 1], pos[3 * j + 2]}; - float d0[3]; - get_mi_vector_dds(d0, pi, pj, box_l, periodic); - - // Loop over every image shift - for (int s = 0; s < d_nShifts; ++s) { - Int3 sh = d_imageShifts[s]; - float dr[3]; - if (sh.x == 0 && sh.y == 0 && sh.z == 0) { - // zero‐shift: use true minimal‐image vector - dr[0] = d0[0]; - dr[1] = d0[1]; - dr[2] = d0[2]; - } else { - // explicit replica: offset the MIC vector - dr[0] = d0[0] + sh.x * box_l[0]; - dr[1] = d0[1] + sh.y * box_l[1]; - dr[2] = d0[2] + sh.z * box_l[2]; - } + float fi[3], ti1[3], ti2[3]; +#ifdef DIPOLE_FIELD_TRACKING + float dfi[3], dfj[3]; +#endif - // Compute interactions for this image - dipole_ia_force(pf, dr, dip + 3 * i, dip + 3 * j, fi, ti, tj); + // --- Minimal‐image branch (no replicas) --- + if (n_replicas == 0) { + for (unsigned int j = i + 1; j < n; ++j) { + // Compute MIC vector + float d0[3]; + get_mi_vector_dds(d0, (float[]){xi, yi, zi}, pos + 3 * j, box_l, + periodic); + // Compute force/torque + dipole_ia_force(pf, d0, mi, dip + 3 * j, fi, ti1, ti2 +#ifdef DIPOLE_FIELD_TRACKING + , + dfi, dfj +#endif + ); - // 1) deposit force/torque on particle j **atomically** - // (many threads i racing to update the same j) + // Deposit on j (atomic) for (int k = 0; k < 3; ++k) { atomicAdd(f + 3 * j + k, -fi[k]); - atomicAdd(torque + 3 * j + k, tj[k]); + atomicAdd(torque + 3 * j + k, ti2[k]); +#ifdef DIPOLE_FIELD_TRACKING + atomicAdd(dip_fld + 3 * j + k, dfj[k]); +#endif } - - // 2) accumulate into thread‐private sums for particle i + // Accumulate on i for (int k = 0; k < 3; ++k) { fsum[k] += fi[k]; - tsum[k] += ti[k]; + tsum[k] += ti1[k]; +#ifdef DIPOLE_FIELD_TRACKING + dfsum[k] += dfi[k]; +#endif + } + } + // Flush i’s result + for (int k = 0; k < 3; ++k) { + atomicAdd(f + 3 * i + k, fsum[k]); + atomicAdd(torque + 3 * i + k, tsum[k]); +#ifdef DIPOLE_FIELD_TRACKING + atomicAdd(dip_fld + 3 * i + k, dfsum[k]); +#endif + } + return; + } + // --- Replica branch (n_replicas > 0) --- + int nx = periodic[0] * n_replicas; + int ny = periodic[1] * n_replicas; + int nz = periodic[2] * n_replicas; + + // Self‐images of i (exclude primary) + for (int dx = -nx; dx <= nx; ++dx) { + for (int dy = -ny; dy <= ny; ++dy) { + for (int dz = -nz; dz <= nz; ++dz) { + if (dx == 0 && dy == 0 && dz == 0) + continue; + float dr[3] = {dx * box_l[0], dy * box_l[1], dz * box_l[2]}; + dipole_ia_force(pf, dr, mi, mi, fi, ti1, ti2 +#ifdef DIPOLE_FIELD_TRACKING + , + dfi, dfj +#endif + ); + for (int k = 0; k < 3; ++k) { + fsum[k] += fi[k]; + tsum[k] += ti1[k]; +#ifdef DIPOLE_FIELD_TRACKING + dfsum[k] += dfi[k]; +#endif + } + } + } + } + + // Pairwise i <-> j over all images + for (unsigned int j = i + 1; j < n; ++j) { + float xj = pos[3 * j + 0], yj = pos[3 * j + 1], zj = pos[3 * j + 2]; + const float *mj = dip + 3 * j; + + for (int dx = -nx; dx <= nx; ++dx) { + for (int dy = -ny; dy <= ny; ++dy) { + for (int dz = -nz; dz <= nz; ++dz) { + float dr[3] = {(xi - xj) + dx * box_l[0], (yi - yj) + dy * box_l[1], + (zi - zj) + dz * box_l[2]}; + dipole_ia_force(pf, dr, mi, mj, fi, ti1, ti2 +#ifdef DIPOLE_FIELD_TRACKING + , + dfi, dfj +#endif + ); + + // Deposit on j (atomic) + for (int k = 0; k < 3; ++k) { + atomicAdd(f + 3 * j + k, -fi[k]); + atomicAdd(torque + 3 * j + k, ti2[k]); +#ifdef DIPOLE_FIELD_TRACKING + atomicAdd(dip_fld + 3 * j + k, dfj[k]); +#endif + fsum[k] += fi[k]; + tsum[k] += ti1[k]; +#ifdef DIPOLE_FIELD_TRACKING + dfsum[k] += dfi[k]; +#endif + } + } } } } @@ -221,6 +289,9 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { atomicAdd(f + 3 * i + k, fsum[k]); atomicAdd(torque + 3 * i + k, tsum[k]); +#ifdef DIPOLE_FIELD_TRACKING + atomicAdd(dip_fld + 3 * i + k, dfsum[k]); +#endif } } // LCOV_EXCL_START @@ -282,44 +353,25 @@ inline void copy_box_data(float **box_l_gpu, int **periodic_gpu, cudaMemcpy(*periodic_gpu, periodic, s_per, cudaMemcpyHostToDevice)); } -// build all shifts with x²+y²+z² ≤ n_replicas² -static std::vector build_shifts(int n_replicas, int periodic[3]) { - std::vector shifts; - int rx = n_replicas * periodic[0]; - int ry = n_replicas * periodic[1]; - int rz = n_replicas * periodic[2]; - int cutoff2 = n_replicas * n_replicas; - shifts.reserve((2 * rx + 1) * (2 * ry + 1) * (2 * rz + 1)); - for (int x = -rx; x <= rx; ++x) - for (int y = -ry; y <= ry; ++y) - for (int z = -rz; z <= rz; ++z) - if (x * x + y * y + z * z <= cutoff2) - shifts.push_back({x, y, z}); - return shifts; -} - void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, - float *dip, float *f, float *torque, + float *dip + +#ifdef DIPOLE_FIELD_TRACKING + , + float *dip_fld +#endif + , + float *f, float *torque, float box_l[3], int periodic[3], int n_replicas) { - unsigned int const bs = 32; + unsigned int const bs = 64; dim3 grid(1, 1, 1); dim3 block(1, 1, 1); if (n == 0) return; - static bool first = true; - if (first) { - auto shifts = build_shifts(n_replicas, periodic); - cudaMemcpyToSymbol(d_imageShifts, shifts.data(), - shifts.size() * sizeof(Int3)); - int nShifts = static_cast(shifts.size()); - cudaMemcpyToSymbol(d_nShifts, &nShifts, sizeof(int)); - first = false; - } - if (n <= bs) { grid.x = 1; block.x = n; @@ -332,8 +384,13 @@ void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, int *periodic_gpu; copy_box_data(&box_l_gpu, &periodic_gpu, box_l, periodic); - KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip, f, - torque, box_l_gpu, periodic_gpu); + KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip +#ifdef DIPOLE_FIELD_TRACKING + , + dip_fld +#endif + , + f, torque, box_l_gpu, periodic_gpu, n_replicas); cudaFree(box_l_gpu); cudaFree(periodic_gpu); } diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh index 8545927328c..079c780a911 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh @@ -32,7 +32,13 @@ void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, float *pos, float *dip, float box_l[3], int periodic[3], float *E); void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, - float *dip, float *f, float *torque, + float *dip +#ifdef DIPOLE_FIELD_TRACKING + , + float *dip_fld +#endif + , + float *f, float *torque, float box_l[3], int periodic[3], int n_replicas); diff --git a/src/core/magnetostatics/dipoles.cpp b/src/core/magnetostatics/dipoles.cpp index d29a1ac7141..df0d1ed8547 100644 --- a/src/core/magnetostatics/dipoles.cpp +++ b/src/core/magnetostatics/dipoles.cpp @@ -181,6 +181,8 @@ struct LongRangeField { actor->dipole_field_at_part(m_particles); } + void operator()(std::shared_ptr const &actor) const {} + template ::value> * = nullptr> void operator()(std::shared_ptr const &) const { diff --git a/src/core/magnetostatics/dipoles.hpp b/src/core/magnetostatics/dipoles.hpp index 3e4d3d96c60..4683d1f7105 100644 --- a/src/core/magnetostatics/dipoles.hpp +++ b/src/core/magnetostatics/dipoles.hpp @@ -72,6 +72,8 @@ using is_solver = std::is_convertible, MagnetostaticsActor>; template struct has_dipole_fields : std::false_type {}; #ifdef DIPOLE_FIELD_TRACKING template <> struct has_dipole_fields : std::true_type {}; +template <> struct has_dipole_fields : std::true_type {}; + #endif // DIPOLE_FIELD_TRACKING } // namespace traits diff --git a/src/core/system/GpuParticleData.cpp b/src/core/system/GpuParticleData.cpp index 75f4df9e26d..fbc2741667a 100644 --- a/src/core/system/GpuParticleData.cpp +++ b/src/core/system/GpuParticleData.cpp @@ -143,6 +143,18 @@ static void add_forces_and_torques(ParticleRange const &particles, i++; } } +#ifdef DIPOLE_FIELD_TRACKING +static void add_dip_fld(ParticleRange const &particles, + std::span dip_fld) { + std::size_t i = 0ul; + for (auto &p : particles) { + for (std::size_t j = 0ul; j < 3ul; j++) { + p.dip_fld()[j] += static_cast(dip_fld[3ul * i + j]); + } + i++; + } +} +#endif /** * @brief Distribute forces to the worker nodes, and add them to the particles. @@ -182,4 +194,26 @@ void GpuParticleData::particles_scatter_forces( } } +#ifdef DIPOLE_FIELD_TRACKING +void GpuParticleData::particles_scatter_dip_fld( + ParticleRange const &particles, std::span host_dip_fld) const { + + auto const size = 3ul * particles.size(); + auto const n_elements = static_cast(size); + + if (::this_node > 0) { + static std::vector buffer_dip_fld; + + buffer_dip_fld.resize(size); + Utils::Mpi::scatter_buffer(buffer_dip_fld.data(), n_elements, ::comm_cart); + add_dip_fld(particles, buffer_dip_fld); + + } else { + Utils::Mpi::scatter_buffer(host_dip_fld.data(), n_elements, ::comm_cart); + + add_dip_fld(particles, host_dip_fld); + } +} +#endif + #endif diff --git a/src/core/system/GpuParticleData.hpp b/src/core/system/GpuParticleData.hpp index 6a933286935..7b6d0209377 100644 --- a/src/core/system/GpuParticleData.hpp +++ b/src/core/system/GpuParticleData.hpp @@ -57,7 +57,8 @@ class GpuParticleData : public System::Leaf { static constexpr std::size_t torque = 2; static constexpr std::size_t q = 3; static constexpr std::size_t dip = 4; - using bitset = std::bitset<5>; + static constexpr std::size_t dip_fld = 5; + using bitset = std::bitset<6>; }; /** @brief Energies that are retrieved from the GPU. */ @@ -101,6 +102,10 @@ class GpuParticleData : public System::Leaf { void particles_scatter_forces(ParticleRange const &particles, std::span host_forces, std::span host_torques) const; +#ifdef DIPOLE_FIELD_TRACKING + void particles_scatter_dip_fld(ParticleRange const &particles, + std::span host_dip_fld) const; +#endif public: GpuParticleData() = default; @@ -115,6 +120,9 @@ class GpuParticleData : public System::Leaf { void enable_property(std::size_t property); void clear_energy_on_device(); void copy_forces_to_host(ParticleRange const &particles, int this_node); +#ifdef DIPOLE_FIELD_TRACKING + void copy_dip_fld_to_host(ParticleRange const &particles, int this_node); +#endif std::size_t n_particles() const; bool has_compatible_device() const; @@ -122,6 +130,9 @@ class GpuParticleData : public System::Leaf { GpuEnergy *get_energy_device() const; float *get_particle_positions_device() const; float *get_particle_forces_device() const; +#ifdef DIPOLE_FIELD_TRACKING + float *get_particle_dip_fld_device() const; +#endif #ifdef ROTATION float *get_particle_torques_device() const; #endif diff --git a/src/core/system/GpuParticleData_cuda.cu b/src/core/system/GpuParticleData_cuda.cu index 1c0be762de4..9aad2b5b31c 100644 --- a/src/core/system/GpuParticleData_cuda.cu +++ b/src/core/system/GpuParticleData_cuda.cu @@ -101,6 +101,10 @@ public: thrust::device_vector particle_data_device; pinned_vector particle_forces_host; thrust::device_vector particle_forces_device; +#ifdef DIPOLE_FIELD_TRACKING + pinned_vector particle_dip_fld_host; + thrust::device_vector particle_dip_fld_device; +#endif #ifdef ROTATION pinned_vector particle_torques_host; thrust::device_vector particle_torques_device; @@ -129,6 +133,15 @@ public: particle_forces_host.begin()); } } +#ifdef DIPOLE_FIELD_TRACKING + void copy_particle_dip_fld_to_host() { + if (not particle_dip_fld_device.empty()) { + thrust::copy(particle_dip_fld_device.begin(), + particle_dip_fld_device.end(), + particle_dip_fld_host.begin()); + } + } +#endif #ifdef ROTATION void copy_particle_torques_to_host() { if (not particle_torques_device.empty()) { @@ -141,6 +154,11 @@ public: std::span get_particle_forces_host_span() { return {particle_forces_host.data(), particle_forces_host.size()}; } +#ifdef DIPOLE_FIELD_TRACKING + std::span get_particle_dip_fld_host_span() { + return {particle_dip_fld_host.data(), particle_dip_fld_host.size()}; + } +#endif #ifdef ROTATION std::span get_particle_torques_host_span() { return {particle_torques_host.data(), particle_torques_host.size()}; @@ -163,6 +181,11 @@ float *GpuParticleData::get_particle_positions_device() const { float *GpuParticleData::get_particle_forces_device() const { return raw_data_pointer(m_data->particle_forces_device); } +#ifdef DIPOLE_FIELD_TRACKING +float *GpuParticleData::get_particle_dip_fld_device() const { + return raw_data_pointer(m_data->particle_dip_fld_device); +} +#endif #ifdef ROTATION float *GpuParticleData::get_particle_torques_device() const { @@ -189,9 +212,16 @@ GpuParticleData::GpuEnergy *GpuParticleData::get_energy_device() const { void GpuParticleData::enable_property(std::size_t property) { m_need_particles_update = true; m_data->m_need[property] = true; +#ifdef DIPOLE_FIELD_TRACKING + if (property != prop::force and property != prop::torque and + property != prop::dip_fld) { + m_split_particle_struct = true; + } +#else if (property != prop::force and property != prop::torque) { m_split_particle_struct = true; } +#endif enable_particle_transfer(); } @@ -222,6 +252,10 @@ void GpuParticleData::Storage::copy_particles_to_device() { resize_or_replace(particle_data_device, n_part); particle_forces_host.resize(3ul * n_part); resize_or_replace(particle_forces_device, 3ul * n_part); +#ifdef DIPOLE_FIELD_TRACKING + particle_dip_fld_host.resize(3ul * n_part); + resize_or_replace(particle_dip_fld_device, 3ul * n_part); +#endif #ifdef ROTATION particle_torques_host.resize(3ul * n_part); resize_or_replace(particle_torques_device, 3ul * n_part); @@ -230,6 +264,10 @@ void GpuParticleData::Storage::copy_particles_to_device() { // zero out device memory for forces and torques cudaMemsetAsync(raw_data_pointer(particle_forces_device), 0x0, byte_size(particle_forces_device), stream[0]); +#ifdef DIPOLE_FIELD_TRACKING + cudaMemsetAsync(raw_data_pointer(particle_dip_fld_device), 0x0, + byte_size(particle_dip_fld_device), stream[0]); +#endif #ifdef ROTATION cudaMemsetAsync(raw_data_pointer(particle_torques_device), 0x0, byte_size(particle_torques_device), stream[0]); @@ -277,6 +315,22 @@ void GpuParticleData::copy_forces_to_host(ParticleRange const &particles, particles_scatter_forces(particles, forces_buffer, torques_buffer); } } +#ifdef DIPOLE_FIELD_TRACKING +void GpuParticleData::copy_dip_fld_to_host(ParticleRange const &particles, + int this_node) { + if (m_communication_enabled) { + // copy results from device memory to host memory + if (this_node == 0) { + m_data->copy_particle_dip_fld_to_host(); + } + + auto dipole_field_buffer = m_data->get_particle_dip_fld_host_span(); + + // add dip_fld to the particles + particles_scatter_dip_fld(particles, dipole_field_buffer); + } +} +#endif void GpuParticleData::clear_energy_on_device() { if (m_communication_enabled) { @@ -426,6 +480,9 @@ void GpuParticleData::Storage::free_device_memory() { }; free_device_vector(particle_data_device); free_device_vector(particle_forces_device); +#ifdef DIPOLE_FIELD_TRACKING + free_device_vector(particle_dip_fld_device); +#endif #ifdef ROTATION free_device_vector(particle_torques_device); #endif diff --git a/testsuite/python/dipole_field_tracking.py b/testsuite/python/dipole_field_tracking.py index aedfb115724..262c6284598 100644 --- a/testsuite/python/dipole_field_tracking.py +++ b/testsuite/python/dipole_field_tracking.py @@ -115,6 +115,14 @@ def test_dds(self): rel_diff = 100. * (time_series[-1] - time_series[0]) / time_series[0] self.assertGreater(np.linalg.norm(rel_diff), 10) + # @utx.skipIfMissingFeatures(["CUDA",]) + def test_dds_gpu(self): + for replicas in [0, 1 ,2]: + solver = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1.,n_replicas=replicas) + self.system.magnetostatics.solver = solver + self.system.integrator.run(steps=1) + for p in self.system.part.all(): + np.testing.assert_allclose(np.copy(p.torque_lab), np.cross(p.dip,p.dip_fld),rtol=1e-9, atol=1e-5) if __name__ == "__main__": ut.main() From 4dde10299ee6917be7a568e5bfb311156dc8f6da Mon Sep 17 00:00:00 2001 From: stekajack Date: Thu, 22 May 2025 13:37:36 +0000 Subject: [PATCH 04/16] cuda guard added --- src/core/magnetostatics/dipoles.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/magnetostatics/dipoles.hpp b/src/core/magnetostatics/dipoles.hpp index 4683d1f7105..5c5dcf86d1a 100644 --- a/src/core/magnetostatics/dipoles.hpp +++ b/src/core/magnetostatics/dipoles.hpp @@ -72,7 +72,9 @@ using is_solver = std::is_convertible, MagnetostaticsActor>; template struct has_dipole_fields : std::false_type {}; #ifdef DIPOLE_FIELD_TRACKING template <> struct has_dipole_fields : std::true_type {}; +#ifdef CUDA template <> struct has_dipole_fields : std::true_type {}; +#endif #endif // DIPOLE_FIELD_TRACKING From 793faafcb1fb7d660716c976c406de235b47a18d Mon Sep 17 00:00:00 2001 From: stekajack Date: Thu, 22 May 2025 15:14:24 +0000 Subject: [PATCH 05/16] remove LongRangeField path --- src/core/energy.cpp | 6 --- .../magnetostatics/dipolar_direct_sum.cpp | 3 ++ src/core/magnetostatics/dipoles.cpp | 37 ++++--------------- src/core/observables/ParticleDipoleFields.hpp | 35 +++++------------- src/python/espressomd/analyze.py | 7 ---- src/script_interface/analysis/Analysis.cpp | 6 --- testsuite/python/dipolar_interface.py | 3 -- testsuite/python/dipole_field_tracking.py | 3 +- testsuite/python/observables.py | 2 +- 9 files changed, 23 insertions(+), 79 deletions(-) diff --git a/src/core/energy.cpp b/src/core/energy.cpp index b1735d14f39..5c50960190b 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -158,10 +158,4 @@ std::optional System::particle_bond_energy(int pid, int bond_id, } } -#ifdef DIPOLE_FIELD_TRACKING -void System::calculate_long_range_fields() { - dipoles.calc_long_range_field(cell_structure->local_particles()); -} -#endif - } // namespace System diff --git a/src/core/magnetostatics/dipolar_direct_sum.cpp b/src/core/magnetostatics/dipolar_direct_sum.cpp index 66fe6ec9319..e894bce3163 100644 --- a/src/core/magnetostatics/dipolar_direct_sum.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum.cpp @@ -368,6 +368,9 @@ void DipolarDirectSum::add_long_range_forces( (*p)->force() += prefactor * fi.f; (*p)->torque() += prefactor * fi.torque; } +#ifdef DIPOLE_FIELD_TRACKING + DipolarDirectSum::dipole_field_at_part(particles); +#endif } /** diff --git a/src/core/magnetostatics/dipoles.cpp b/src/core/magnetostatics/dipoles.cpp index df0d1ed8547..e71df758e96 100644 --- a/src/core/magnetostatics/dipoles.cpp +++ b/src/core/magnetostatics/dipoles.cpp @@ -138,6 +138,14 @@ struct LongRangeForce { actor->add_long_range_forces(); } #endif +#ifdef DIPOLE_FIELD_TRACKING + template ::value> * = nullptr> + void operator()(std::shared_ptr const &) const { + runtimeErrorMsg() << "Dipoles field calculation not implemented by " + << "dipolar method " << Utils::demangle(); + } +#endif }; struct LongRangeEnergy { @@ -171,27 +179,6 @@ struct LongRangeEnergy { #endif }; -#ifdef DIPOLE_FIELD_TRACKING -struct LongRangeField { - ParticleRange const &m_particles; - explicit LongRangeField(ParticleRange const &particles) - : m_particles(particles) {} - - void operator()(std::shared_ptr const &actor) const { - actor->dipole_field_at_part(m_particles); - } - - void operator()(std::shared_ptr const &actor) const {} - - template ::value> * = nullptr> - void operator()(std::shared_ptr const &) const { - runtimeErrorMsg() << "Dipoles field calculation not implemented by " - << "dipolar method " << Utils::demangle(); - } -}; -#endif - void Solver::calc_pressure_long_range() const { if (impl->solver) { runtimeWarningMsg() << "pressure calculated, but pressure not implemented."; @@ -211,13 +198,5 @@ double Solver::calc_energy_long_range(ParticleRange const &particles) const { return 0.; } -#ifdef DIPOLE_FIELD_TRACKING -void Solver::calc_long_range_field(ParticleRange const &particles) const { - if (impl->solver) { - std::visit(LongRangeField(particles), *impl->solver); - } -} -#endif - } // namespace Dipoles #endif // DIPOLES diff --git a/src/core/observables/ParticleDipoleFields.hpp b/src/core/observables/ParticleDipoleFields.hpp index 131bbd12777..fd19b56385a 100644 --- a/src/core/observables/ParticleDipoleFields.hpp +++ b/src/core/observables/ParticleDipoleFields.hpp @@ -17,36 +17,21 @@ * along with this program. If not, see . */ -#pragma once - -#include "config/config.hpp" +#ifndef OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP +#define OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP #include "PidObservable.hpp" -#include "system/System.hpp" - -#include namespace Observables { -/** Extract particle dipole fields. - * For \f$n\f$ particles, return \f$3 n\f$ dipole fields ordered as - * \f$(h_d1_x, h_d1_y, h_d1_z, \dots, h_dn_x, h_dn_y, h_dn_z)\f$. +/** Extract dip_flds. + * For \f$n\f$ particles, return \f$3 n\f$ OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP + * ordered as + * \f$(dip_fld_x^1, dip_fld_y^1, dip_fld_z^1, \dots, dip_fld_x^n, dip_fld_y^n, + * dip_fld_z^n)\f$. */ -class ParticleDipoleFields - : public ParticleObservable { -public: - using ParticleObservable< - ParticleObservables::DipoleFields>::ParticleObservable; - std::vector - evaluate(boost::mpi::communicator const &comm, - ParticleReferenceRange const &local_particles, - const ParticleObservables::traits &traits) const override { -#ifdef DIPOLE_FIELD_TRACKING - System::get_system().calculate_long_range_fields(); -#endif - return ParticleObservable::evaluate( - comm, local_particles, traits); - } -}; +using ParticleDipoleFields = + ParticleObservable; } // namespace Observables +#endif \ No newline at end of file diff --git a/src/python/espressomd/analyze.py b/src/python/espressomd/analyze.py index c3af629f464..7820ef96202 100644 --- a/src/python/espressomd/analyze.py +++ b/src/python/espressomd/analyze.py @@ -541,13 +541,6 @@ def particle_bond_energy(self, particle, bond): return self.call_method("particle_bond_energy", pid=particle.id, bond_id=interaction._bond_id, partners=partners) - def dipole_fields(self): - """ - Calculate the total dipole field on each particle. - """ - assert_features("DIPOLE_FIELD_TRACKING") - self.call_method("calc_long_range_fields") - def dpd_stress(self): assert_features("DPD") return np.reshape(self.call_method("dpd_stress"), (3, 3)) diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index 17c4f1ad498..08953a7f029 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -135,12 +135,6 @@ Variant Analysis::do_call_method(std::string const &name, auto const local = system.particle_bond_energy(pid, bond_id, partners); return Utils::Mpi::reduce_optional(context()->get_comm(), local); } -#ifdef DIPOLE_FIELD_TRACKING - if (name == "calc_long_range_fields") { - get_system().calculate_long_range_fields(); - return {}; - } -#endif if (name == "particle_neighbor_pids") { auto &system = get_system(); system.on_observable_calc(); diff --git a/testsuite/python/dipolar_interface.py b/testsuite/python/dipolar_interface.py index 1f6500f0bb5..1a17de6ced5 100644 --- a/testsuite/python/dipolar_interface.py +++ b/testsuite/python/dipolar_interface.py @@ -115,9 +115,6 @@ def test_exceptions_non_p3m(self): ["DIPOLAR_DIRECT_SUM", "DIPOLE_FIELD_TRACKING"]) and has_gpu: ddsg = DDSG(prefactor=1.) self.system.magnetostatics.solver = ddsg - with self.assertRaisesRegex(Exception, "Dipoles field calculation not implemented by dipolar method DipolarDirectSumGpu"): - self.system.part.add(pos=(0.2, 0.2, 0.2), dip=(0.0, 0.0, 1.0)) - self.system.analysis.dipole_fields() self.system.part.clear() self.system.magnetostatics.clear() # check it's safe to resize the box, i.e. there are no currently diff --git a/testsuite/python/dipole_field_tracking.py b/testsuite/python/dipole_field_tracking.py index 262c6284598..2e3e7b37b61 100644 --- a/testsuite/python/dipole_field_tracking.py +++ b/testsuite/python/dipole_field_tracking.py @@ -98,7 +98,6 @@ def test_dds(self): solver = espressomd.magnetostatics.DipolarDirectSumCpu(prefactor=1.) self.system.magnetostatics.solver = solver self.system.integrator.run(steps=0) - self.system.analysis.dipole_fields() slice_data = [(x.id, x.pos, x.dip) for x in self.system.part.all()] dip_fields_obs = espressomd.observables.ParticleDipoleFields( ids=self.system.part.all().id) @@ -117,7 +116,7 @@ def test_dds(self): # @utx.skipIfMissingFeatures(["CUDA",]) def test_dds_gpu(self): - for replicas in [0, 1 ,2]: + for replicas in [0, 1]: solver = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1.,n_replicas=replicas) self.system.magnetostatics.solver = solver self.system.integrator.run(steps=1) diff --git a/testsuite/python/observables.py b/testsuite/python/observables.py index 1385c7bfd43..6395d2a77a2 100644 --- a/testsuite/python/observables.py +++ b/testsuite/python/observables.py @@ -53,7 +53,7 @@ class Observables(ut.TestCase): partcls.dip = np.random.random((N_PART, 3)) - .3 if espressomd.has_features(["DIPOLE_FIELD_TRACKING"]): - system.analysis.dipole_fields() + partcls.dip_fld = np.random.random((N_PART, 3)) - .3 if espressomd.has_features(["ROTATION"]): partcls.omega_body = np.random.random((N_PART, 3)) - .5 From feb112e7c3d73a1cae853b7374059c37b38cf430 Mon Sep 17 00:00:00 2001 From: stekajack Date: Fri, 23 May 2025 09:31:04 +0000 Subject: [PATCH 06/16] add missing feature checks in tests --- testsuite/python/dipolar_direct_summation.py | 5 ++++- testsuite/python/dipole_field_tracking.py | 8 +++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/testsuite/python/dipolar_direct_summation.py b/testsuite/python/dipolar_direct_summation.py index cceec180717..2ef70e2440e 100644 --- a/testsuite/python/dipolar_direct_summation.py +++ b/testsuite/python/dipolar_direct_summation.py @@ -296,7 +296,10 @@ def test_inner_loop_consistency_cpu(self): np.testing.assert_allclose(node_01_energy, node_00_energy, **tol) np.testing.assert_allclose(node_01_forces, node_00_forces, **tol) np.testing.assert_allclose(node_01_torques, node_00_torques, **tol) - + + @utx.skipIfMissingGPU() + @ut.skipIf(system.cell_system.get_state()["n_nodes"] == 1, + "only runs for 2 or more MPI ranks") def test_inner_loop_consistency_gpu(self): system = self.system system.periodicity = [True, True, True] diff --git a/testsuite/python/dipole_field_tracking.py b/testsuite/python/dipole_field_tracking.py index 2e3e7b37b61..f84bc4de05c 100644 --- a/testsuite/python/dipole_field_tracking.py +++ b/testsuite/python/dipole_field_tracking.py @@ -114,14 +114,16 @@ def test_dds(self): rel_diff = 100. * (time_series[-1] - time_series[0]) / time_series[0] self.assertGreater(np.linalg.norm(rel_diff), 10) - # @utx.skipIfMissingFeatures(["CUDA",]) + @utx.skipIfMissingGPU() def test_dds_gpu(self): for replicas in [0, 1]: - solver = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1.,n_replicas=replicas) + solver = espressomd.magnetostatics.DipolarDirectSumGpu( + prefactor=1., n_replicas=replicas) self.system.magnetostatics.solver = solver self.system.integrator.run(steps=1) for p in self.system.part.all(): - np.testing.assert_allclose(np.copy(p.torque_lab), np.cross(p.dip,p.dip_fld),rtol=1e-9, atol=1e-5) + np.testing.assert_allclose(np.copy(p.torque_lab), np.cross( + p.dip,p.dip_fld), rtol=1e-9, atol=1e-5) if __name__ == "__main__": ut.main() From 99ea83dc5a71509898c9eefbe0e4895a34e78048 Mon Sep 17 00:00:00 2001 From: stekajack Date: Fri, 23 May 2025 12:13:14 +0000 Subject: [PATCH 07/16] formatting --- testsuite/python/dipole_field_tracking.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/python/dipole_field_tracking.py b/testsuite/python/dipole_field_tracking.py index f84bc4de05c..3960a205d11 100644 --- a/testsuite/python/dipole_field_tracking.py +++ b/testsuite/python/dipole_field_tracking.py @@ -123,7 +123,7 @@ def test_dds_gpu(self): self.system.integrator.run(steps=1) for p in self.system.part.all(): np.testing.assert_allclose(np.copy(p.torque_lab), np.cross( - p.dip,p.dip_fld), rtol=1e-9, atol=1e-5) + p.dip, p.dip_fld), rtol=1e-9, atol=1e-5) if __name__ == "__main__": ut.main() From 74d1e91b16c628d97f8d1bd42dff5a4df9c05e72 Mon Sep 17 00:00:00 2001 From: stekajack Date: Fri, 23 May 2025 12:17:56 +0000 Subject: [PATCH 08/16] formatting2 --- testsuite/python/dipole_field_tracking.py | 1 + 1 file changed, 1 insertion(+) diff --git a/testsuite/python/dipole_field_tracking.py b/testsuite/python/dipole_field_tracking.py index 3960a205d11..730d7f869ed 100644 --- a/testsuite/python/dipole_field_tracking.py +++ b/testsuite/python/dipole_field_tracking.py @@ -125,5 +125,6 @@ def test_dds_gpu(self): np.testing.assert_allclose(np.copy(p.torque_lab), np.cross( p.dip, p.dip_fld), rtol=1e-9, atol=1e-5) + if __name__ == "__main__": ut.main() From c7cd617042ade7558193b6a84305c88feb41d8d4 Mon Sep 17 00:00:00 2001 From: stekajack Date: Tue, 10 Jun 2025 15:19:30 +0000 Subject: [PATCH 09/16] addresed implicit conversion warnign and fixed caliper test fail --- .../dipolar_direct_sum_gpu_cuda.cu | 9 ++++++--- testsuite/python/caliper.py | 17 +++++++++++++++-- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu index dab708f3059..d8f4237fba6 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu @@ -231,7 +231,9 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int dz = -nz; dz <= nz; ++dz) { if (dx == 0 && dy == 0 && dz == 0) continue; - float dr[3] = {dx * box_l[0], dy * box_l[1], dz * box_l[2]}; + float dr[3] = {static_cast(dx) * box_l[0], + static_cast(dy) * box_l[1], + static_cast(dz) * box_l[2]}; dipole_ia_force(pf, dr, mi, mi, fi, ti1, ti2 #ifdef DIPOLE_FIELD_TRACKING , @@ -257,8 +259,9 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int dx = -nx; dx <= nx; ++dx) { for (int dy = -ny; dy <= ny; ++dy) { for (int dz = -nz; dz <= nz; ++dz) { - float dr[3] = {(xi - xj) + dx * box_l[0], (yi - yj) + dy * box_l[1], - (zi - zj) + dz * box_l[2]}; + float dr[3] = {(xi - xj) + static_cast(dx) * box_l[0], + (yi - yj) + static_cast(dy) * box_l[1], + (zi - zj) + static_cast(dz) * box_l[2]}; dipole_ia_force(pf, dr, mi, mj, fi, ti1, ti2 #ifdef DIPOLE_FIELD_TRACKING , diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py index 67aedbb11bc..09ffd881250 100644 --- a/testsuite/python/caliper.py +++ b/testsuite/python/caliper.py @@ -31,6 +31,7 @@ calculate_forces copy_particles_to_GPU init_forces + invalidate_dip_fld calc_long_range_forces short_range_loop copy_forces_from_GPU @@ -38,6 +39,7 @@ calculate_forces copy_particles_to_GPU init_forces + invalidate_dip_fld calc_long_range_forces short_range_loop copy_forces_from_GPU @@ -52,6 +54,8 @@ class Test(ut.TestCase): @utx.skipIfMissingFeatures(["P3M", "WCA"]) def test_runtime_report(self): has_cuda = espressomd.has_features(["CUDA"]) + has_dipfld = espressomd.has_features(["DIPOLE_FIELD_TRACKING"]) + script = str(pathlib.Path(__file__).parent / "caliper_child.py") my_env = os.environ.copy() my_env["CALI_CONFIG"] = "runtime-report" @@ -69,8 +73,17 @@ def test_runtime_report(self): self.assertEqual(lines[0].split(), header.split(), msg=f"Caliper summary should start with '{header}'") labels = [line[:30].strip() for line in lines[1:]] - labels_ref = [x.strip() for x in EXPECTED_LABELS.strip().split("\n") - if "GPU" not in x.upper() or has_cuda] + # build expected labels, skipping GPU-only and dip_fld if not enabled + labels_ref = [] + for line in EXPECTED_LABELS.strip().split("\n"): + label = line.strip() + # skip GPU-only entries if CUDA not present + if "GPU" in label.upper() and not has_cuda: + continue + # skip our dip_fld entry if the feature wasn't compiled in + if label == "invalidate_dip_fld" and not has_dipfld: + continue + labels_ref.append(label) self.assertEqual(labels[:len(labels_ref)], labels_ref, msg=f"Caliper returned this summary:\n{stderr}") From 76d6c76d0e8bfbfc75638657833bf07da38436cf Mon Sep 17 00:00:00 2001 From: stekajack Date: Thu, 25 Sep 2025 15:57:46 +0000 Subject: [PATCH 10/16] update preprocessor directives --- .../magnetostatics/dipolar_direct_sum.cpp | 2 +- .../magnetostatics/dipolar_direct_sum_gpu.cpp | 6 +-- .../dipolar_direct_sum_gpu_cuda.cu | 46 ++++++++++++------- .../dipolar_direct_sum_gpu_cuda.cuh | 2 +- src/core/magnetostatics/dipoles.hpp | 1 - src/core/system/GpuParticleData.cpp | 4 +- src/core/virtual_sites.cpp | 3 +- 7 files changed, 38 insertions(+), 26 deletions(-) diff --git a/src/core/magnetostatics/dipolar_direct_sum.cpp b/src/core/magnetostatics/dipolar_direct_sum.cpp index 70913b570f5..3e57d46b126 100644 --- a/src/core/magnetostatics/dipolar_direct_sum.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum.cpp @@ -368,7 +368,7 @@ void DipolarDirectSum::add_long_range_forces( (*p)->force() += prefactor * fi.f; (*p)->torque() += prefactor * fi.torque; } -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING DipolarDirectSum::dipole_field_at_part(particles); #endif } diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index 639378d6076..9ddb018cfbb 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -51,7 +51,7 @@ void DipolarDirectSumGpu::on_activation() const { gpu_particle_data.enable_property(GpuParticleData::prop::torque); gpu_particle_data.enable_property(GpuParticleData::prop::pos); gpu_particle_data.enable_property(GpuParticleData::prop::dip); -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING gpu_particle_data.enable_property(GpuParticleData::prop::dip_fld); #endif } @@ -71,12 +71,12 @@ void DipolarDirectSumGpu::add_long_range_forces() const { auto const torques_device = gpu.get_particle_torques_device(); auto const positions_device = gpu.get_particle_positions_device(); auto const dipoles_device = gpu.get_particle_dipoles_device(); -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING auto const dipole_fields_device = gpu.get_particle_dip_fld_device(); #endif DipolarDirectSum_kernel_wrapper_force( static_cast(prefactor), npart, positions_device, dipoles_device -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , dipole_fields_device #endif diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu index 3d8ede3b97f..846f5fd0487 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu @@ -59,13 +59,25 @@ __device__ inline void get_mi_vector_dds(float res[3], float const a[3], __device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1, float const *dip2, float *f1, float *torque1, float *torque2 -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , float *dip_fld1, float *dip_fld2 #endif ) { // Powers of distance auto const r_sq = scalar_product(dr, dr); + if (r_sq == 0.0f) { + f1[0] = f1[1] = f1[2] = 0.0f; +#ifdef ROTATION + torque1[0] = torque1[1] = torque1[2] = 0.0f; + torque2[0] = torque2[1] = torque2[2] = 0.0f; +#endif +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING + dip_fld1[0] = dip_fld1[1] = dip_fld1[2] = 0.0f; + dip_fld2[0] = dip_fld2[1] = dip_fld2[2] = 0.0f; +#endif + return; + } auto const r_sq_inv = 1.0f / r_sq; auto const r_inv = rsqrtf(r_sq); auto const r3_inv = 1.0f / r_sq * r_inv; @@ -78,7 +90,7 @@ __device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1, auto const pe3 = scalar_product(dip2, dr); auto const pe4 = 3.0f * r5_inv; -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING auto const rep1 = pe3 * pe4; auto const rep2 = pe2 * pe4; dip_fld1[0] = pf * (rep1 * dr[0] - dip2[0] * r3_inv); @@ -148,7 +160,7 @@ __device__ float dipole_ia_energy(float pf, float const *r1, float const *r2, __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, const float *pos, const float *dip -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , float *dip_fld #endif @@ -169,12 +181,12 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, // Per‐thread accumulators float fsum[3] = {0.0f, 0.0f, 0.0f}; float tsum[3] = {0.0f, 0.0f, 0.0f}; -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING float dfsum[3] = {0.0f, 0.0f, 0.0f}; #endif float fi[3], ti1[3], ti2[3]; -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING float dfi[3], dfj[3]; #endif @@ -187,7 +199,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, periodic); // Compute force/torque dipole_ia_force(pf, d0, mi, dip + 3 * j, fi, ti1, ti2 -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , dfi, dfj #endif @@ -197,7 +209,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { atomicAdd(f + 3 * j + k, -fi[k]); atomicAdd(torque + 3 * j + k, ti2[k]); -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING atomicAdd(dip_fld + 3 * j + k, dfj[k]); #endif } @@ -205,7 +217,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { fsum[k] += fi[k]; tsum[k] += ti1[k]; -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING dfsum[k] += dfi[k]; #endif } @@ -214,7 +226,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { atomicAdd(f + 3 * i + k, fsum[k]); atomicAdd(torque + 3 * i + k, tsum[k]); -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING atomicAdd(dip_fld + 3 * i + k, dfsum[k]); #endif } @@ -235,7 +247,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, static_cast(dy) * box_l[1], static_cast(dz) * box_l[2]}; dipole_ia_force(pf, dr, mi, mi, fi, ti1, ti2 -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , dfi, dfj #endif @@ -243,7 +255,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { fsum[k] += fi[k]; tsum[k] += ti1[k]; -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING dfsum[k] += dfi[k]; #endif } @@ -263,7 +275,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, (yi - yj) + static_cast(dy) * box_l[1], (zi - zj) + static_cast(dz) * box_l[2]}; dipole_ia_force(pf, dr, mi, mj, fi, ti1, ti2 -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , dfi, dfj #endif @@ -273,12 +285,12 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { atomicAdd(f + 3 * j + k, -fi[k]); atomicAdd(torque + 3 * j + k, ti2[k]); -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING atomicAdd(dip_fld + 3 * j + k, dfj[k]); #endif fsum[k] += fi[k]; tsum[k] += ti1[k]; -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING dfsum[k] += dfi[k]; #endif } @@ -292,7 +304,7 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, for (int k = 0; k < 3; ++k) { atomicAdd(f + 3 * i + k, fsum[k]); atomicAdd(torque + 3 * i + k, tsum[k]); -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING atomicAdd(dip_fld + 3 * i + k, dfsum[k]); #endif } @@ -359,7 +371,7 @@ inline void copy_box_data(float **box_l_gpu, int **periodic_gpu, void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, float *dip -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , float *dip_fld #endif @@ -388,7 +400,7 @@ void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, copy_box_data(&box_l_gpu, &periodic_gpu, box_l, periodic); KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , dip_fld #endif diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh index cefd2d0ad2e..9a212df3963 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh @@ -33,7 +33,7 @@ void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, float *pos, int periodic[3], float *E); void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, float *dip -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING , float *dip_fld #endif diff --git a/src/core/magnetostatics/dipoles.hpp b/src/core/magnetostatics/dipoles.hpp index dabb89e4c08..07806b4109e 100644 --- a/src/core/magnetostatics/dipoles.hpp +++ b/src/core/magnetostatics/dipoles.hpp @@ -77,7 +77,6 @@ template <> struct has_dipole_fields : std::true_type {}; #endif #endif // ESPRESSO_DIPOLE_FIELD_TRACKING - } // namespace traits } // namespace Dipoles diff --git a/src/core/system/GpuParticleData.cpp b/src/core/system/GpuParticleData.cpp index 3a918670894..fac4dc59a02 100644 --- a/src/core/system/GpuParticleData.cpp +++ b/src/core/system/GpuParticleData.cpp @@ -143,7 +143,7 @@ static void add_forces_and_torques(ParticleRange const &particles, i++; } } -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING static void add_dip_fld(ParticleRange const &particles, std::span dip_fld) { std::size_t i = 0ul; @@ -194,7 +194,7 @@ void GpuParticleData::particles_scatter_forces( } } -#ifdef DIPOLE_FIELD_TRACKING +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING void GpuParticleData::particles_scatter_dip_fld( ParticleRange const &particles, std::span host_dip_fld) const { diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index 0b4db1ac9c9..ae6bc941f10 100644 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -107,8 +107,9 @@ calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, constexpr auto const *function_name = location.function_name(); constexpr auto *error_msg = "%s: component %u: %f instead of %f\n"; Utils::Quaternion qtemp = relate_to_quat * quat; + auto const epsilon = 1e-10; for (unsigned int i = 0; i < 4; i++) { - if (fabs(qtemp[i] - quat_director[i]) != 0.) { + if (fabs(qtemp[i] - quat_director[i]) >= epsilon) { fprintf(stderr, error_msg, function_name, i, qtemp[i], quat_director[i]); } } From bf01249e12f2492f287ec67c7516c02b3a736c7a Mon Sep 17 00:00:00 2001 From: stekajack Date: Fri, 26 Sep 2025 13:06:04 +0000 Subject: [PATCH 11/16] fix caliper --- src/core/forces.cpp | 2 +- testsuite/python/caliper.py | 24 ++++++++++++------------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 7a4e2d87465..db9818947ee 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -137,7 +137,7 @@ void init_forces_and_thermostat(System::System const &system) { } #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING void invalidate_dip_fld(const CellStructure &cell_structure) { -#ifdef CALIPER +#ifdef ESPRESSO_CALIPER CALI_CXX_MARK_FUNCTION; #endif diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py index 1ad9c5583cc..7a57e21528d 100644 --- a/testsuite/python/caliper.py +++ b/testsuite/python/caliper.py @@ -32,17 +32,17 @@ Initial Force Calculation calculate_forces copy_particles_to_GPU - invalidate_dip_fld {'convert particles AoS to SoA' if HAS_CABANA else ''} - init_forces_and_thermost + invalidate_dip_fld + init_forces_and_thermostat calc_long_range_forces {'parallel short range' if HAS_CABANA else 'short_range_loop'} copy_forces_from_GPU Integration loop calculate_forces copy_particles_to_GPU - invalidate_dip_fld {'convert particles AoS to SoA' if HAS_CABANA else ''} + invalidate_dip_fld init_forces_and_thermostat calc_long_range_forces {'parallel short range' if HAS_CABANA else 'short_range_loop'} @@ -77,18 +77,18 @@ def test_runtime_report(self): self.assertEqual(lines[0].split(), header.split(), msg=f"Caliper summary should start with '{header}'") - labels = [line[:30].strip() for line in lines[1:]] - # build expected labels, skipping GPU-only and dip_fld if not enabled + labels = [line[:36].rstrip() for line in lines[1:]] + labels_ref = [] - for line in EXPECTED_LABELS.strip().split("\n"): - label = line.strip() - # skip GPU-only entries if CUDA not present - if "GPU" in label.upper() and not has_cuda: + for x in EXPECTED_LABELS.strip().split("\n"): + x = x.rstrip() + if not x: + continue + if "GPU" in x.upper() and not has_cuda: continue - # skip our dip_fld entry if the feature wasn't compiled in - if label == "invalidate_dip_fld" and not has_dipfld: + if x.strip() == "invalidate_dip_fld" and not has_dipfld: continue - labels_ref.append(label) + labels_ref.append(x) self.assertEqual(labels[:len(labels_ref)], labels_ref, msg=f"Caliper returned this summary:\n{stderr}") From 991fdd634096337627219fe7da2b5668f0ed906e Mon Sep 17 00:00:00 2001 From: stekajack Date: Fri, 26 Sep 2025 13:39:42 +0000 Subject: [PATCH 12/16] fix: reorder 'convert particles AoS to SoA' in caliper.py --- testsuite/python/caliper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py index 7a57e21528d..a315ff5d64c 100644 --- a/testsuite/python/caliper.py +++ b/testsuite/python/caliper.py @@ -32,8 +32,8 @@ Initial Force Calculation calculate_forces copy_particles_to_GPU - {'convert particles AoS to SoA' if HAS_CABANA else ''} invalidate_dip_fld + {'convert particles AoS to SoA' if HAS_CABANA else ''} init_forces_and_thermostat calc_long_range_forces {'parallel short range' if HAS_CABANA else 'short_range_loop'} @@ -41,8 +41,8 @@ Integration loop calculate_forces copy_particles_to_GPU - {'convert particles AoS to SoA' if HAS_CABANA else ''} invalidate_dip_fld + {'convert particles AoS to SoA' if HAS_CABANA else ''} init_forces_and_thermostat calc_long_range_forces {'parallel short range' if HAS_CABANA else 'short_range_loop'} From ce164825c0a83855f51e7e76bb952fdb7ec590a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 9 Oct 2025 15:46:19 +0200 Subject: [PATCH 13/16] Patching --- doc/sphinx/magnetostatics.rst | 4 ++-- src/core/forces.cpp | 4 ---- src/core/observables/ParticleDipoleFields.hpp | 13 +++++-------- src/core/system/GpuParticleData.hpp | 2 -- src/core/system/GpuParticleData_cuda.cu | 12 ++++-------- src/core/virtual_sites.cpp | 2 +- src/python/espressomd/magnetostatics.py | 8 ++++++-- testsuite/python/caliper.py | 19 ++----------------- 8 files changed, 20 insertions(+), 44 deletions(-) diff --git a/doc/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst index e2b3e923988..7c5decd1cd7 100644 --- a/doc/sphinx/magnetostatics.rst +++ b/doc/sphinx/magnetostatics.rst @@ -156,10 +156,10 @@ attach it to the system. The only required parameter is the prefactor dds = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1) system.magnetostatics.solver = dds -The CPU implementation has an optional argument ``n_replicas`` which +Both implementations have an optional argument ``n_replicas`` which adds periodic copies to the system along periodic directions. In that case, the minimum image convention is no longer used. -Additionally, enabling the ``DIPOLE_FIELDS_TRACKING`` feature enables the CPU +Additionally, enabling the ``DIPOLE_FIELD_TRACKING`` feature enables either implementation to calculate the total dipole field at the position of each magnetic particle in the primary simulation box. These values are stored in the particle handle's ``dip_fld`` property and can be accessed directly or diff --git a/src/core/forces.cpp b/src/core/forces.cpp index db9818947ee..091ba8ea162 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -137,10 +137,6 @@ void init_forces_and_thermostat(System::System const &system) { } #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING void invalidate_dip_fld(const CellStructure &cell_structure) { -#ifdef ESPRESSO_CALIPER - CALI_CXX_MARK_FUNCTION; -#endif - cell_structure.for_each_local_particle( [](Particle &p) { p.dip_fld() = {0., 0., 0.}; }); } diff --git a/src/core/observables/ParticleDipoleFields.hpp b/src/core/observables/ParticleDipoleFields.hpp index fd19b56385a..eb5fd09fd6c 100644 --- a/src/core/observables/ParticleDipoleFields.hpp +++ b/src/core/observables/ParticleDipoleFields.hpp @@ -17,21 +17,18 @@ * along with this program. If not, see . */ -#ifndef OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP -#define OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP +#pragma once #include "PidObservable.hpp" namespace Observables { -/** Extract dip_flds. - * For \f$n\f$ particles, return \f$3 n\f$ OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP - * ordered as - * \f$(dip_fld_x^1, dip_fld_y^1, dip_fld_z^1, \dots, dip_fld_x^n, dip_fld_y^n, - * dip_fld_z^n)\f$. +/** + * @brief Extract dip_flds. + * For \f$n\f$ particles, return \f$3 n\f$ dipole fields ordered as + * \f$(h_d1_x, h_d1_y, h_d1_z, \dots, h_dn_x, h_dn_y, h_dn_z)\f$. */ using ParticleDipoleFields = ParticleObservable; } // namespace Observables -#endif \ No newline at end of file diff --git a/src/core/system/GpuParticleData.hpp b/src/core/system/GpuParticleData.hpp index 28fb86bf279..16bf8c44c26 100644 --- a/src/core/system/GpuParticleData.hpp +++ b/src/core/system/GpuParticleData.hpp @@ -130,12 +130,10 @@ class GpuParticleData : public System::Leaf { GpuEnergy *get_energy_device() const; float *get_particle_positions_device() const; float *get_particle_forces_device() const; - #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING float *get_particle_dip_fld_device() const; #endif #ifdef ESPRESSO_ROTATION - float *get_particle_torques_device() const; #endif #ifdef ESPRESSO_DIPOLES diff --git a/src/core/system/GpuParticleData_cuda.cu b/src/core/system/GpuParticleData_cuda.cu index 89f65866d23..6af89ab3922 100644 --- a/src/core/system/GpuParticleData_cuda.cu +++ b/src/core/system/GpuParticleData_cuda.cu @@ -101,13 +101,11 @@ public: thrust::device_vector particle_data_device; pinned_vector particle_forces_host; thrust::device_vector particle_forces_device; - #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING pinned_vector particle_dip_fld_host; thrust::device_vector particle_dip_fld_device; #endif #ifdef ESPRESSO_ROTATION - pinned_vector particle_torques_host; thrust::device_vector particle_torques_device; #endif @@ -145,8 +143,8 @@ public: } } #endif -#ifdef ESPRESSO_ROTATION +#ifdef ESPRESSO_ROTATION void copy_particle_torques_to_host() { if (not particle_torques_device.empty()) { thrust::copy(particle_torques_device.begin(), @@ -164,8 +162,8 @@ public: return {particle_dip_fld_host.data(), particle_dip_fld_host.size()}; } #endif -#ifdef ESPRESSO_ROTATION +#ifdef ESPRESSO_ROTATION std::span get_particle_torques_host_span() { return {particle_torques_host.data(), particle_torques_host.size()}; } @@ -263,8 +261,8 @@ void GpuParticleData::Storage::copy_particles_to_device() { particle_dip_fld_host.resize(3ul * n_part); resize_or_replace(particle_dip_fld_device, 3ul * n_part); #endif -#ifdef ESPRESSO_ROTATION +#ifdef ESPRESSO_ROTATION particle_torques_host.resize(3ul * n_part); resize_or_replace(particle_torques_device, 3ul * n_part); #endif @@ -277,8 +275,8 @@ void GpuParticleData::Storage::copy_particles_to_device() { cudaMemsetAsync(raw_data_pointer(particle_dip_fld_device), 0x0, byte_size(particle_dip_fld_device), stream[0]); #endif -#ifdef ESPRESSO_ROTATION +#ifdef ESPRESSO_ROTATION cudaMemsetAsync(raw_data_pointer(particle_torques_device), 0x0, byte_size(particle_torques_device), stream[0]); #endif @@ -490,12 +488,10 @@ void GpuParticleData::Storage::free_device_memory() { }; free_device_vector(particle_data_device); free_device_vector(particle_forces_device); - #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING free_device_vector(particle_dip_fld_device); #endif #ifdef ESPRESSO_ROTATION - free_device_vector(particle_torques_device); #endif free_device_pointer(particle_pos_device); diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index ae6bc941f10..da27bb4fb6e 100644 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -107,7 +107,7 @@ calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, constexpr auto const *function_name = location.function_name(); constexpr auto *error_msg = "%s: component %u: %f instead of %f\n"; Utils::Quaternion qtemp = relate_to_quat * quat; - auto const epsilon = 1e-10; + constexpr auto epsilon = 1e-10; for (unsigned int i = 0; i < 4; i++) { if (fabs(qtemp[i] - quat_director[i]) >= epsilon) { fprintf(stderr, error_msg, function_name, i, qtemp[i], quat_director[i]); diff --git a/src/python/espressomd/magnetostatics.py b/src/python/espressomd/magnetostatics.py index e41776cbac0..b939a95c0c2 100644 --- a/src/python/espressomd/magnetostatics.py +++ b/src/python/espressomd/magnetostatics.py @@ -235,8 +235,10 @@ class DipolarDirectSumGpu(MagnetostaticInteraction): Calculate magnetostatic interactions by direct summation over all pairs. See :ref:`Dipolar direct sum` for more details. - If the system has periodic boundaries, the minimum image convention - is applied in the respective directions. + If the system has periodic boundaries, the minimum image convention is + applied in the respective directions when no replicas are used. When + replicas are used, ``n_replicas`` copies of the system are taken into + account in the respective directions, and a spherical cutoff is applied. This is the GPU version of :class:`espressomd.magnetostatics.DipolarDirectSumCpu` but uses floating point precision. @@ -248,6 +250,8 @@ class DipolarDirectSumGpu(MagnetostaticInteraction): ---------- prefactor : :obj:`float` Magnetostatics prefactor (:math:`\\mu_0/(4\\pi)`) + n_replicas : :obj:`int` + Number of replicas to be taken into account at periodic boundaries. """ _so_name = "Dipoles::DipolarDirectSumGpu" diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py index a315ff5d64c..30771e10a02 100644 --- a/testsuite/python/caliper.py +++ b/testsuite/python/caliper.py @@ -32,7 +32,6 @@ Initial Force Calculation calculate_forces copy_particles_to_GPU - invalidate_dip_fld {'convert particles AoS to SoA' if HAS_CABANA else ''} init_forces_and_thermostat calc_long_range_forces @@ -41,7 +40,6 @@ Integration loop calculate_forces copy_particles_to_GPU - invalidate_dip_fld {'convert particles AoS to SoA' if HAS_CABANA else ''} init_forces_and_thermostat calc_long_range_forces @@ -58,8 +56,6 @@ class Test(ut.TestCase): @utx.skipIfMissingFeatures(["P3M", "WCA"]) def test_runtime_report(self): has_cuda = espressomd.has_features(["CUDA"]) - has_dipfld = espressomd.has_features(["DIPOLE_FIELD_TRACKING"]) - script = str(pathlib.Path(__file__).parent / "caliper_child.py") my_env = os.environ.copy() my_env["CALI_CONFIG"] = "runtime-report" @@ -76,20 +72,9 @@ def test_runtime_report(self): header = "Path\tMin time/rank\tMax time/rank\tAvg time/rank\tTime %" self.assertEqual(lines[0].split(), header.split(), msg=f"Caliper summary should start with '{header}'") - labels = [line[:36].rstrip() for line in lines[1:]] - - labels_ref = [] - for x in EXPECTED_LABELS.strip().split("\n"): - x = x.rstrip() - if not x: - continue - if "GPU" in x.upper() and not has_cuda: - continue - if x.strip() == "invalidate_dip_fld" and not has_dipfld: - continue - labels_ref.append(x) - + labels_ref = [x.rstrip() for x in EXPECTED_LABELS.strip().split("\n") + if x.rstrip() and ("GPU" not in x.upper() or has_cuda)] self.assertEqual(labels[:len(labels_ref)], labels_ref, msg=f"Caliper returned this summary:\n{stderr}") From 181abe9caee419b3a24b1126fccdcc74d520ac96 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 10 Oct 2025 17:21:35 +0200 Subject: [PATCH 14/16] Reduce code duplication --- src/core/forces.cpp | 15 +- .../magnetostatics/dipolar_direct_sum_gpu.cpp | 12 +- .../dipolar_direct_sum_gpu_cuda.cu | 156 +++++++----------- .../dipolar_direct_sum_gpu_cuda.cuh | 29 ++-- src/core/magnetostatics/dipoles.cpp | 8 - src/core/magnetostatics/dipoles.hpp | 9 - src/script_interface/analysis/Analysis.cpp | 1 - testsuite/python/CMakeLists.txt | 2 +- testsuite/python/dipolar_direct_summation.py | 47 ++---- testsuite/python/dipolar_interface.py | 6 - testsuite/python/dipole_field_tracking.py | 6 +- 11 files changed, 103 insertions(+), 188 deletions(-) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 091ba8ea162..e5ee9b8c521 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -135,12 +135,6 @@ void init_forces_and_thermostat(System::System const &system) { // Initialize ghost forces (unchanged) init_forces_ghosts(cell_structure); } -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING -void invalidate_dip_fld(const CellStructure &cell_structure) { - cell_structure.for_each_local_particle( - [](Particle &p) { p.dip_fld() = {0., 0., 0.}; }); -} -#endif void init_forces_ghosts(const CellStructure &cell_structure) { cell_structure.for_each_ghost_particle( @@ -160,6 +154,13 @@ static void force_capping(CellStructure &cell_structure, double force_cap) { } } +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING +static void reinit_dip_fld(CellStructure const &cell_structure) { + cell_structure.for_each_local_particle( + [](Particle &p) { p.dip_fld() = {0., 0., 0.}; }); +} +#endif + void System::System::calculate_forces() { #ifdef ESPRESSO_CALIPER CALI_CXX_MARK_FUNCTION; @@ -190,7 +191,7 @@ void System::System::calculate_forces() { #endif #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING // reset dipole field - invalidate_dip_fld(*cell_structure); + reinit_dip_fld(*cell_structure); #endif // Use combined function instead of two separate calls diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index 9ddb018cfbb..6381fffa696 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -73,15 +73,13 @@ void DipolarDirectSumGpu::add_long_range_forces() const { auto const dipoles_device = gpu.get_particle_dipoles_device(); #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING auto const dipole_fields_device = gpu.get_particle_dip_fld_device(); +#else + float const *dipole_fields_device{nullptr}; #endif DipolarDirectSum_kernel_wrapper_force( - static_cast(prefactor), npart, positions_device, dipoles_device -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - dipole_fields_device -#endif - , - forces_device, torques_device, box, periodicity, n_replicas); + static_cast(prefactor), npart, positions_device, dipoles_device, + dipole_fields_device, forces_device, torques_device, box, periodicity, + n_replicas); } void DipolarDirectSumGpu::long_range_energy() const { diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu index 846f5fd0487..a5c6ac6c5e8 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu @@ -56,14 +56,12 @@ __device__ inline void get_mi_vector_dds(float res[3], float const a[3], } } -__device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1, - float const *dip2, float *f1, float *torque1, - float *torque2 -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - float *dip_fld1, float *dip_fld2 -#endif -) { +__device__ void dipole_ia_force(float const pf, float const dr[3], + float const *const dip1, + float const *const dip2, float *const f1, + float *const torque1, float *const torque2, + [[maybe_unused]] float *const dip_fld1, + [[maybe_unused]] float *const dip_fld2) { // Powers of distance auto const r_sq = scalar_product(dr, dr); if (r_sq == 0.0f) { @@ -158,25 +156,20 @@ __device__ float dipole_ia_energy(float pf, float const *r1, float const *r2, } // LCOV_EXCL_STOP -__global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, - const float *pos, const float *dip -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - float *dip_fld -#endif - , - float *f, float *torque, - const float box_l[3], - const int periodic[3], - int n_replicas) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; +__global__ void +DipolarDirectSum_kernel_force(float const pf, unsigned int n, + float const *const pos, float const *const dip, + [[maybe_unused]] float *dip_fld, float *const f, + float *const torque, float const box_l[3], + int const periodic[3], int const n_replicas) { + unsigned int const i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= n) return; // Load particle i - float xi = pos[3 * i + 0], yi = pos[3 * i + 1], zi = pos[3 * i + 2]; + auto const xi{pos[3u * i + 0u]}, yi{pos[3u * i + 1u]}, zi{pos[3u * i + 2u]}; - const float *mi = dip + 3 * i; + auto const *const mi = dip + 3u * i; // Per‐thread accumulators float fsum[3] = {0.0f, 0.0f, 0.0f}; @@ -188,33 +181,29 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, float fi[3], ti1[3], ti2[3]; #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING float dfi[3], dfj[3]; +#else + float *const dfi{nullptr}, *const dfj{nullptr}; #endif // --- Minimal‐image branch (no replicas) --- if (n_replicas == 0) { - for (unsigned int j = i + 1; j < n; ++j) { + for (unsigned int j = i + 1u; j < n; ++j) { // Compute MIC vector float d0[3]; - get_mi_vector_dds(d0, (float[]){xi, yi, zi}, pos + 3 * j, box_l, - periodic); + get_mi_vector_dds(d0, pos + 3u * i, pos + 3u * j, box_l, periodic); // Compute force/torque - dipole_ia_force(pf, d0, mi, dip + 3 * j, fi, ti1, ti2 -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - dfi, dfj -#endif - ); + dipole_ia_force(pf, d0, mi, dip + 3u * j, fi, ti1, ti2, dfi, dfj); // Deposit on j (atomic) - for (int k = 0; k < 3; ++k) { - atomicAdd(f + 3 * j + k, -fi[k]); - atomicAdd(torque + 3 * j + k, ti2[k]); + for (unsigned int k = 0u; k < 3u; ++k) { + atomicAdd(f + 3u * j + k, -fi[k]); + atomicAdd(torque + 3u * j + k, ti2[k]); #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - atomicAdd(dip_fld + 3 * j + k, dfj[k]); + atomicAdd(dip_fld + 3u * j + k, dfj[k]); #endif } // Accumulate on i - for (int k = 0; k < 3; ++k) { + for (unsigned int k = 0u; k < 3u; ++k) { fsum[k] += fi[k]; tsum[k] += ti1[k]; #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING @@ -222,20 +211,20 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, #endif } } - // Flush i’s result - for (int k = 0; k < 3; ++k) { - atomicAdd(f + 3 * i + k, fsum[k]); - atomicAdd(torque + 3 * i + k, tsum[k]); + // Flush i's result + for (unsigned int k = 0u; k < 3u; ++k) { + atomicAdd(f + 3u * i + k, fsum[k]); + atomicAdd(torque + 3u * i + k, tsum[k]); #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - atomicAdd(dip_fld + 3 * i + k, dfsum[k]); + atomicAdd(dip_fld + 3u * i + k, dfsum[k]); #endif } return; } // --- Replica branch (n_replicas > 0) --- - int nx = periodic[0] * n_replicas; - int ny = periodic[1] * n_replicas; - int nz = periodic[2] * n_replicas; + int const nx = periodic[0] * n_replicas; + int const ny = periodic[1] * n_replicas; + int const nz = periodic[2] * n_replicas; // Self‐images of i (exclude primary) for (int dx = -nx; dx <= nx; ++dx) { @@ -246,13 +235,8 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, float dr[3] = {static_cast(dx) * box_l[0], static_cast(dy) * box_l[1], static_cast(dz) * box_l[2]}; - dipole_ia_force(pf, dr, mi, mi, fi, ti1, ti2 -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - dfi, dfj -#endif - ); - for (int k = 0; k < 3; ++k) { + dipole_ia_force(pf, dr, mi, mi, fi, ti1, ti2, dfi, dfj); + for (unsigned int k = 0u; k < 3u; ++k) { fsum[k] += fi[k]; tsum[k] += ti1[k]; #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING @@ -264,9 +248,9 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, } // Pairwise i <-> j over all images - for (unsigned int j = i + 1; j < n; ++j) { - float xj = pos[3 * j + 0], yj = pos[3 * j + 1], zj = pos[3 * j + 2]; - const float *mj = dip + 3 * j; + for (unsigned int j = i + 1u; j < n; ++j) { + auto const xj{pos[3u * j + 0u]}, yj{pos[3u * j + 1u]}, zj{pos[3u * j + 2u]}; + auto const *mj = dip + 3u * j; for (int dx = -nx; dx <= nx; ++dx) { for (int dy = -ny; dy <= ny; ++dy) { @@ -274,19 +258,14 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, float dr[3] = {(xi - xj) + static_cast(dx) * box_l[0], (yi - yj) + static_cast(dy) * box_l[1], (zi - zj) + static_cast(dz) * box_l[2]}; - dipole_ia_force(pf, dr, mi, mj, fi, ti1, ti2 -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - dfi, dfj -#endif - ); + dipole_ia_force(pf, dr, mi, mj, fi, ti1, ti2, dfi, dfj); // Deposit on j (atomic) - for (int k = 0; k < 3; ++k) { - atomicAdd(f + 3 * j + k, -fi[k]); - atomicAdd(torque + 3 * j + k, ti2[k]); + for (unsigned int k = 0u; k < 3u; ++k) { + atomicAdd(f + 3u * j + k, -fi[k]); + atomicAdd(torque + 3u * j + k, ti2[k]); #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - atomicAdd(dip_fld + 3 * j + k, dfj[k]); + atomicAdd(dip_fld + 3u * j + k, dfj[k]); #endif fsum[k] += fi[k]; tsum[k] += ti1[k]; @@ -299,13 +278,13 @@ __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n, } } - // Add i’s total to global memory **atomically**, + // Add i's total to global memory **atomically**, // in case any other asynchronous update might race - for (int k = 0; k < 3; ++k) { - atomicAdd(f + 3 * i + k, fsum[k]); - atomicAdd(torque + 3 * i + k, tsum[k]); + for (unsigned int k = 0u; k < 3u; ++k) { + atomicAdd(f + 3u * i + k, fsum[k]); + atomicAdd(torque + 3u * i + k, tsum[k]); #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - atomicAdd(dip_fld + 3 * i + k, dfsum[k]); + atomicAdd(dip_fld + 3u * i + k, dfsum[k]); #endif } } @@ -327,7 +306,8 @@ __device__ void dds_sumReduction(float *input, float *sum) { // LCOV_EXCL_STOP __global__ void DipolarDirectSum_kernel_energy(float pf, unsigned int n, - float *pos, float *dip, + float const *const pos, + float const *const dip, float box_l[3], int periodic[3], float *energySum) { @@ -368,17 +348,12 @@ inline void copy_box_data(float **box_l_gpu, int **periodic_gpu, cudaMemcpy(*periodic_gpu, periodic, s_per, cudaMemcpyHostToDevice)); } -void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, - float *dip - -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - float *dip_fld -#endif - , - float *f, float *torque, - float box_l[3], int periodic[3], - int n_replicas) { +void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, + float const *const pos, + float const *const dip, + float *dip_fld, float *f, + float *torque, float box_l[3], + int periodic[3], int n_replicas) { unsigned int const bs = 64; dim3 grid(1, 1, 1); @@ -399,20 +374,17 @@ void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, int *periodic_gpu; copy_box_data(&box_l_gpu, &periodic_gpu, box_l, periodic); - KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - dip_fld -#endif - , - f, torque, box_l_gpu, periodic_gpu, n_replicas); + KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip, + dip_fld, f, torque, box_l_gpu, periodic_gpu, n_replicas); cudaFree(box_l_gpu); cudaFree(periodic_gpu); } -void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, float *pos, - float *dip, float box_l[3], - int periodic[3], float *E) { +void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, + float const *const pos, + float const *const dip, + float box_l[3], int periodic[3], + float *E) { unsigned int const bs = 512; dim3 grid(1, 1, 1); @@ -453,4 +425,4 @@ void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, float *pos, cuda_safe_mem(cudaFree(periodic_gpu)); } -#endif +#endif // ESPRESSO_DIPOLAR_DIRECT_SUM diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh index 9a212df3963..4e3d1b1b178 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cuh @@ -23,23 +23,16 @@ #ifdef ESPRESSO_DIPOLAR_DIRECT_SUM -// simple 3-int struct for image shifts -struct Int3 { - int x, y, z; -}; - -void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, float *pos, - float *dip, float box_l[3], - int periodic[3], float *E); -void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos, - float *dip -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - , - float *dip_fld -#endif - , - float *f, float *torque, - float box_l[3], int periodic[3], - int n_replicas); +void DipolarDirectSum_kernel_wrapper_energy(float k, unsigned int n, + float const *const pos, + float const *const dip, + float box_l[3], int periodic[3], + float *E); +void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, + float const *const pos, + float const *const dip, + float *dip_fld, float *f, + float *torque, float box_l[3], + int periodic[3], int n_replicas); #endif // ESPRESSO_DIPOLAR_DIRECT_SUM diff --git a/src/core/magnetostatics/dipoles.cpp b/src/core/magnetostatics/dipoles.cpp index 3bd2a148c63..9f08195cb7a 100644 --- a/src/core/magnetostatics/dipoles.cpp +++ b/src/core/magnetostatics/dipoles.cpp @@ -138,14 +138,6 @@ struct LongRangeForce { actor->add_long_range_forces(); } #endif -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - template ::value> * = nullptr> - void operator()(std::shared_ptr const &) const { - runtimeErrorMsg() << "Dipoles field calculation not implemented by " - << "dipolar method " << Utils::demangle(); - } -#endif }; struct LongRangeEnergy { diff --git a/src/core/magnetostatics/dipoles.hpp b/src/core/magnetostatics/dipoles.hpp index 07806b4109e..ca7cddfdc70 100644 --- a/src/core/magnetostatics/dipoles.hpp +++ b/src/core/magnetostatics/dipoles.hpp @@ -68,15 +68,6 @@ namespace traits { template using is_solver = std::is_convertible, MagnetostaticsActor>; -/** @brief The dipolar method supports dipole fields calculation. */ -template struct has_dipole_fields : std::false_type {}; -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING -template <> struct has_dipole_fields : std::true_type {}; -#ifdef ESPRESSO_CUDA -template <> struct has_dipole_fields : std::true_type {}; -#endif -#endif // ESPRESSO_DIPOLE_FIELD_TRACKING - } // namespace traits } // namespace Dipoles diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index 949d3740d06..a0d83942450 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -135,7 +135,6 @@ Variant Analysis::do_call_method(std::string const &name, auto const local = system.particle_bond_energy(pid, bond_id, partners); return Utils::Mpi::reduce_optional(context()->get_comm(), local); } - if (name == "particle_neighbor_pids") { auto &system = get_system(); system.on_observable_calc(); diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index b63d5776237..075a66ad438 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -318,7 +318,7 @@ python_test(FILE coulomb_cloud_wall_duplicated.py MAX_NUM_PROC 4 GPU_SLOTS 3) python_test(FILE collision_detection.py MAX_NUM_PROC 4) python_test(FILE collision_detection_interface.py MAX_NUM_PROC 2) python_test(FILE lj.py MAX_NUM_PROC 4) -python_test(FILE dipole_field_tracking.py MAX_NUM_PROC 2) +python_test(FILE dipole_field_tracking.py MAX_NUM_PROC 2 GPU_SLOTS 1) python_test(FILE pairs.py MAX_NUM_PROC 4) python_test(FILE polymer_linear.py MAX_NUM_PROC 4) python_test(FILE polymer_diamond.py MAX_NUM_PROC 4) diff --git a/testsuite/python/dipolar_direct_summation.py b/testsuite/python/dipolar_direct_summation.py index 2ef70e2440e..73bd9a8009c 100644 --- a/testsuite/python/dipolar_direct_summation.py +++ b/testsuite/python/dipolar_direct_summation.py @@ -262,57 +262,30 @@ def test_min_image_convention_gpu(self): @ut.skipIf(system.cell_system.get_state()["n_nodes"] == 1, "only runs for 2 or more MPI ranks") def test_inner_loop_consistency_cpu(self): - system = self.system - system.periodicity = [True, True, True] - tol = {"atol": 1e-10, "rtol": 1e-10} - p1 = system.part.add(pos=[0., 0., 0.], dip=[0., 0., 1.], - rotation=[True, True, True]) - p2 = system.part.add(pos=[1., 0., 0.], dip=[0., 0., 1.], - rotation=[True, True, True]) - for n_replicas in [0, 1]: - system.magnetostatics.clear() - solver = espressomd.magnetostatics.DipolarDirectSumCpu( - prefactor=1., n_replicas=n_replicas) - system.magnetostatics.solver = solver - - # intra-node calculation - p1.pos = [system.box_l[0] / 2. - 0.1, 0., 2.] - p2.pos = [system.box_l[0] / 2. + 0.1, 0., 0.] - system.integrator.run(steps=0, recalc_forces=True) - assert p1.node != p2.node - node_01_energy = system.analysis.energy()["dipolar"] - node_01_forces = np.copy(system.part.all().f) - node_01_torques = np.copy(system.part.all().torque_lab) - - # inter-node calculation - p1.pos = [0.1, 0., 2.] - p2.pos = [0.3, 0., 0.] - system.integrator.run(steps=0, recalc_forces=True) - assert p1.node == p2.node - node_00_energy = system.analysis.energy()["dipolar"] - node_00_forces = np.copy(system.part.all().f) - node_00_torques = np.copy(system.part.all().torque_lab) - - np.testing.assert_allclose(node_01_energy, node_00_energy, **tol) - np.testing.assert_allclose(node_01_forces, node_00_forces, **tol) - np.testing.assert_allclose(node_01_torques, node_00_torques, **tol) + self.check_inner_loop_consistency( + solver=espressomd.magnetostatics.DipolarDirectSumCpu, + tol={"atol": 1e-10, "rtol": 1e-10}) + @utx.skipIfMissingFeatures("DIPOLAR_DIRECT_SUM") @utx.skipIfMissingGPU() @ut.skipIf(system.cell_system.get_state()["n_nodes"] == 1, "only runs for 2 or more MPI ranks") def test_inner_loop_consistency_gpu(self): + self.check_inner_loop_consistency( + solver=espressomd.magnetostatics.DipolarDirectSumGpu, + tol={"atol": 1e-6, "rtol": 1e-6}) + + def check_inner_loop_consistency(self, solver, tol): system = self.system system.periodicity = [True, True, True] - tol = {"atol": 1e-6, "rtol": 1e-6} p1 = system.part.add(pos=[0., 0., 0.], dip=[0., 0., 1.], rotation=[True, True, True]) p2 = system.part.add(pos=[1., 0., 0.], dip=[0., 0., 1.], rotation=[True, True, True]) for n_replicas in [0, 1]: system.magnetostatics.clear() - solver = espressomd.magnetostatics.DipolarDirectSumGpu( + system.magnetostatics.solver = solver( prefactor=1., n_replicas=n_replicas) - system.magnetostatics.solver = solver # intra-node calculation p1.pos = [system.box_l[0] / 2. - 0.1, 0., 2.] diff --git a/testsuite/python/dipolar_interface.py b/testsuite/python/dipolar_interface.py index 1a17de6ced5..00a7f1bb140 100644 --- a/testsuite/python/dipolar_interface.py +++ b/testsuite/python/dipolar_interface.py @@ -111,12 +111,6 @@ def test_exceptions_non_p3m(self): mdlc = MDLC(gap_size=1., maxPWerror=1e-5, actor=ddsr) self.system.magnetostatics.solver = mdlc self.assertIsNone(self.system.magnetostatics.solver) - if espressomd.has_features( - ["DIPOLAR_DIRECT_SUM", "DIPOLE_FIELD_TRACKING"]) and has_gpu: - ddsg = DDSG(prefactor=1.) - self.system.magnetostatics.solver = ddsg - self.system.part.clear() - self.system.magnetostatics.clear() # check it's safe to resize the box, i.e. there are no currently # active sanity check in the core self.system.change_volume_and_rescale_particles(10., "y") diff --git a/testsuite/python/dipole_field_tracking.py b/testsuite/python/dipole_field_tracking.py index 730d7f869ed..83f258a1117 100644 --- a/testsuite/python/dipole_field_tracking.py +++ b/testsuite/python/dipole_field_tracking.py @@ -115,6 +115,7 @@ def test_dds(self): self.assertGreater(np.linalg.norm(rel_diff), 10) @utx.skipIfMissingGPU() + @utx.skipIfMissingFeatures(["DIPOLAR_DIRECT_SUM"]) def test_dds_gpu(self): for replicas in [0, 1]: solver = espressomd.magnetostatics.DipolarDirectSumGpu( @@ -122,8 +123,9 @@ def test_dds_gpu(self): self.system.magnetostatics.solver = solver self.system.integrator.run(steps=1) for p in self.system.part.all(): - np.testing.assert_allclose(np.copy(p.torque_lab), np.cross( - p.dip, p.dip_fld), rtol=1e-9, atol=1e-5) + np.testing.assert_allclose( + np.copy(p.torque_lab), np.cross(p.dip, p.dip_fld), + rtol=1e-9, atol=1e-5) if __name__ == "__main__": From 94d502ca3afca572fc8ffde46c3de32ae70cbf76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 10 Oct 2025 17:45:50 +0200 Subject: [PATCH 15/16] Fix regression --- src/core/magnetostatics/dipolar_direct_sum_gpu.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index 6381fffa696..3f611e936ba 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -74,7 +74,7 @@ void DipolarDirectSumGpu::add_long_range_forces() const { #ifdef ESPRESSO_DIPOLE_FIELD_TRACKING auto const dipole_fields_device = gpu.get_particle_dip_fld_device(); #else - float const *dipole_fields_device{nullptr}; + float *dipole_fields_device{nullptr}; #endif DipolarDirectSum_kernel_wrapper_force( static_cast(prefactor), npart, positions_device, dipoles_device, From e5e68c71c8c9abda391413fe52110948a9816660 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 10 Oct 2025 18:03:50 +0200 Subject: [PATCH 16/16] fixup --- src/core/magnetostatics/dipolar_direct_sum_gpu.cpp | 7 +++---- testsuite/python/dipolar_interface.py | 2 ++ 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index 3f611e936ba..a57b8ed30c3 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -71,10 +71,9 @@ void DipolarDirectSumGpu::add_long_range_forces() const { auto const torques_device = gpu.get_particle_torques_device(); auto const positions_device = gpu.get_particle_positions_device(); auto const dipoles_device = gpu.get_particle_dipoles_device(); -#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING - auto const dipole_fields_device = gpu.get_particle_dip_fld_device(); -#else float *dipole_fields_device{nullptr}; +#ifdef ESPRESSO_DIPOLE_FIELD_TRACKING + dipole_fields_device = gpu.get_particle_dip_fld_device(); #endif DipolarDirectSum_kernel_wrapper_force( static_cast(prefactor), npart, positions_device, dipoles_device, @@ -101,4 +100,4 @@ void DipolarDirectSumGpu::long_range_energy() const { periodicity, energy_device); } -#endif +#endif // ESPRESSO_DIPOLAR_DIRECT_SUM diff --git a/testsuite/python/dipolar_interface.py b/testsuite/python/dipolar_interface.py index 00a7f1bb140..b02ccfa530a 100644 --- a/testsuite/python/dipolar_interface.py +++ b/testsuite/python/dipolar_interface.py @@ -90,6 +90,8 @@ def test_exceptions_non_p3m(self): ddsg = DDSG(prefactor=1.) with self.assertRaisesRegex(ValueError, "Parameter 'actor' of type Dipoles::DipolarDirectSumGpu isn't supported by DLC"): MDLC(gap_size=2., maxPWerror=0.1, actor=ddsg) + with self.assertRaisesRegex(ValueError, "Parameter 'n_replicas' must be >= 0"): + DDSG(prefactor=1., n_replicas=-2) with self.assertRaisesRegex(RuntimeError, "Parameter 'actor' is missing"): MDLC(gap_size=2., maxPWerror=0.1) with self.assertRaisesRegex(RuntimeError, "Parameter 'n_replica' is not a valid parameter"):