From c867159dba7b4e653061250fe791f92c454c7fc4 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Tue, 14 Nov 2023 03:02:19 +0100
Subject: [PATCH] Add half-neighbor lists variant
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
examples/dem.py | 2 +-
src/pairs/ir/apply.py | 5 +++--
src/pairs/sim/interaction.py | 26 +++++++++++++++++++++++---
src/pairs/transformations/loops.py | 5 ++++-
4 files changed, 31 insertions(+), 7 deletions(-)
diff --git a/examples/dem.py b/examples/dem.py
index ecd4889..60f86c0 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -137,7 +137,6 @@ psim = pairs.simulation(
double_prec=True,
use_contact_history=True)
-#psim = pairs.simulation("dem", debug=True, timesteps=timeSteps)
psim.add_position('position')
psim.add_property('mass', pairs.real(), 1.0)
psim.add_property('linear_velocity', pairs.vector())
@@ -177,6 +176,7 @@ psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
'pi': math.pi,
'infinity': math.inf })
+#psim.compute_half()
psim.build_neighbor_lists(linkedCellWidth + skin)
psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
diff --git a/src/pairs/ir/apply.py b/src/pairs/ir/apply.py
index e2658e7..43c21e4 100644
--- a/src/pairs/ir/apply.py
+++ b/src/pairs/ir/apply.py
@@ -7,6 +7,7 @@ from pairs.ir.properties import Property
from pairs.ir.scalars import ScalarOp
from pairs.ir.vectors import Vector
from pairs.ir.types import Types
+from pairs.sim.flags import Flags
from pairs.sim.lowerable import FinalLowerable, Lowerable
@@ -45,8 +46,8 @@ class Apply(Lowerable):
def lower(self):
if self.sim._compute_half:
for _ in Filter(self.sim,
- ScalarOp.and_op(j < self.sim.nlocal,
- ScalarOp.cmp(self.sim.particle_flags[j] & Flags.Fixed, 0))):
+ 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)
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index 3fcf52f..fb5becc 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -59,6 +59,7 @@ class NeighborFor:
ncells = self.cell_lists.ncells
particle_cell = self.cell_lists.particle_cell
cell_particles = self.cell_lists.cell_particles
+ particle_shape = self.sim.particle_shape
nshapes = self.cell_lists.nshapes
for shape in range(self.sim.max_shapes()):
@@ -67,17 +68,36 @@ class NeighborFor:
for _ in Filter(self.sim, ScalarOp.and_op(neigh_cell > 0, neigh_cell < ncells)):
start = sum([nshapes[neigh_cell][s] for s in range(shape)], 0)
-
for cell_particle in For(self.sim, start, start + nshapes[neigh_cell][shape]):
particle_id = cell_particles[neigh_cell][cell_particle]
- for _ in Filter(self.sim, ScalarOp.neq(particle_id, self.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))
+
+ else:
+ condition = ScalarOp.neq(particle_id, self.particle)
+
+ for _ in Filter(self.sim, condition):
yield Neighbor(self.sim, cell_particle, neigh_cell, particle_id, shape)
# Infinite particles
start = sum([nshapes[0][s] for s in range(shape)], 0)
for inf_id in For(self.sim, start, start + nshapes[0][shape]):
inf_particle = cell_particles[0][inf_id]
- for _ in Filter(self.sim, ScalarOp.neq(inf_particle, self.particle)):
+
+ 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))
+
+ else:
+ condition = ScalarOp.neq(inf_particle, self.particle)
+
+ for _ in Filter(self.sim, condition):
yield Neighbor(self.sim, inf_id, 0, inf_particle, shape)
else:
diff --git a/src/pairs/transformations/loops.py b/src/pairs/transformations/loops.py
index 829f6f7..238effd 100644
--- a/src/pairs/transformations/loops.py
+++ b/src/pairs/transformations/loops.py
@@ -42,7 +42,10 @@ class LICM(Mutator):
Matrix,
MatrixOp,
PropertyAccess,
- Quaternion,
+ # TODO: This is commented to avoid errors from the CUDA Explicit Euler kernel used
+ # for the DEM test cases. Allow lifting of Vector, Quaternion and Matrix until the
+ # kernel loops only (maybe other elements as well)
+ #Quaternion,
QuaternionOp,
ScalarOp,
Select,
--
GitLab