From d6f0b279f9da4549b8edaaafad05b6a6d633fe63 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Thu, 16 Nov 2023 00:11:58 +0100
Subject: [PATCH] Fix half-neighbor lists
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
examples/lj.py | 3 ++-
src/pairs/ir/apply.py | 4 ++--
src/pairs/sim/interaction.py | 18 +++++++++---------
3 files changed, 13 insertions(+), 12 deletions(-)
diff --git a/examples/lj.py b/examples/lj.py
index 97b66a4..214966e 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -5,7 +5,7 @@ import sys
def lj(i, j):
sr2 = 1.0 / squared_distance(i, j)
sr6 = sr2 * sr2 * sr2 * sigma6[i, j]
- apply(force, delta(i, j) * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j])
+ apply(force, delta(i, j) * (48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j]))
def euler(i):
linear_velocity[i] += dt * force[i] / mass[i]
@@ -36,6 +36,7 @@ psim.add_feature_property('type', 'sigma6', pairs.real(), [epsilon for i in rang
psim.set_domain([0.0, 0.0, 0.0, 53.747078, 53.747078, 53.747078])
psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'linear_velocity', 'flags'], pairs.point_mass())
psim.reneighbor_every(20)
+#psim.compute_half()
psim.build_neighbor_lists(cutoff_radius + skin)
psim.vtk_output(f"output/lj_{target}")
psim.compute(lj, cutoff_radius)
diff --git a/src/pairs/ir/apply.py b/src/pairs/ir/apply.py
index 43c21e4..6ae09dd 100644
--- a/src/pairs/ir/apply.py
+++ b/src/pairs/ir/apply.py
@@ -44,11 +44,11 @@ class Apply(Lowerable):
@pairs_inline
def lower(self):
+ Assign(self.sim, self._reduction_variable, self._reduction_variable + self._expr)
+
if self.sim._compute_half:
for _ in Filter(self.sim,
ScalarOp.and_op(self._j < self.sim.nlocal,
ScalarOp.cmp(self.sim.particle_flags[self._j] & Flags.Fixed, 0))):
Assign(self.sim, self._prop[self._j], self._prop[self._j] - self._expr)
-
- Assign(self.sim, self._reduction_variable, self._reduction_variable + self._expr)
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index fb5becc..ff2f5c2 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -72,10 +72,8 @@ class NeighborFor:
particle_id = cell_particles[neigh_cell][cell_particle]
if self.sim._compute_half:
- shape_cond = particle_shape[particle_id] > shape
- condition = ScalarOp.or_op(shape_cond,
- ScalarOp.and_op(ScalarOp.not_op(shape_cond),
- self.particle < particle_id))
+ shape_cond = particle_shape[particle_id] > self.sim.get_shape_id(shape)
+ condition = ScalarOp.or_op(shape_cond, self.particle < particle_id)
else:
condition = ScalarOp.neq(particle_id, self.particle)
@@ -89,10 +87,8 @@ class NeighborFor:
inf_particle = cell_particles[0][inf_id]
if self.sim._compute_half:
- shape_cond = particle_shape[inf_particle] > shape
- condition = ScalarOp.or_op(shape_cond,
- ScalarOp.and_op(ScalarOp.not_op(shape_cond),
- self.particle < inf_particle))
+ shape_cond = particle_shape[inf_particle] > self.sim.get_shape_id(shape)
+ condition = ScalarOp.or_op(shape_cond, self.particle < inf_particle)
else:
condition = ScalarOp.neq(inf_particle, self.particle)
@@ -266,7 +262,11 @@ class ParticleInteraction(Lowerable):
prop_reductions[prop] = prop_reductions[prop] + reduction
for prop, reduction in prop_reductions.items():
- Assign(self.sim, prop[i], reduction)
+ if self.sim._compute_half:
+ Assign(self.sim, prop[i], prop[i] + reduction)
+
+ else:
+ Assign(self.sim, prop[i], reduction)
else:
raise Exception("Interactions among more than two particles are currently not supported.")
--
GitLab