From acba635515e571c38f0f14c5eb15d4a7d914d649 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Sat, 18 Nov 2023 02:19:06 +0100
Subject: [PATCH] Implement contact history for linked cells

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
 Makefile                    |  9 +++++----
 examples/dem.py             |  1 +
 src/pairs/mapping/funcs.py  | 34 ++++++++++++++++++++++++++++++++--
 src/pairs/sim/simulation.py |  4 +++-
 4 files changed, 41 insertions(+), 7 deletions(-)

diff --git a/Makefile b/Makefile
index 343a2dc..e083715 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 .PHONY: all build clean
 
-TESTCASE=lj
+TESTCASE=dem
 PYCMD=python3
 CC=mpicxx
 NVCC=nvcc
@@ -12,7 +12,8 @@ CPU_SRC="$(TESTCASE).cpp"
 CPU_BIN="$(TESTCASE)_cpu"
 GPU_SRC="$(TESTCASE).cu"
 GPU_BIN="$(TESTCASE)_gpu"
-DEBUG_FLAGS="-DDEBUG"
+DEBUG_FLAGS=
+#DEBUG_FLAGS="-DDEBUG"
 
 all: clean build $(CPU_BIN) $(GPU_BIN)
 	@echo "Everything was done!"
@@ -42,10 +43,10 @@ $(OBJ_PATH)/dummy.o: runtime/devices/dummy.cpp
 
 # Targets
 $(CPU_BIN): $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o
-	$(CC) -O3 -o $(CPU_BIN) $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o -DDEBUG
+	$(CC) -O3 -o $(CPU_BIN) $(CPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o $(OBJ_PATH)/dummy.o ${DEBUG_FLAGS}
 
 $(GPU_BIN): $(GPU_SRC) $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o 
-	$(NVCC) -c -o $(OBJ_PATH)/cuda_runtime.o runtime/devices/cuda.cu -DDEBUG
+	$(NVCC) -c -o $(OBJ_PATH)/cuda_runtime.o runtime/devices/cuda.cu ${DEBUG_FLAGS}
 	$(NVCC) -c -o $(OBJ_PATH)/$(GPU_BIN).o $(GPU_SRC) -DDEBUG
 	$(CC) -o $(GPU_BIN) $(OBJ_PATH)/$(GPU_BIN).o $(OBJ_PATH)/cuda_runtime.o $(OBJ_PATH)/pairs.o $(OBJ_PATH)/regular_6d_stencil.o -lcudart -L$(CUDA_PATH)/lib64
 
diff --git a/examples/dem.py b/examples/dem.py
index b8630d2..1bb5ea9 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -177,6 +177,7 @@ psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
                                      'infinity': math.inf })
 
 #psim.compute_half()
+#psim.build_cell_lists(linkedCellWidth)
 psim.build_neighbor_lists(linkedCellWidth + skin)
 psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
 
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index d27f86e..0d93008 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -3,9 +3,10 @@ import inspect
 from pairs.ir.assign import Assign
 from pairs.ir.branches import Branch, Filter
 from pairs.ir.lit import Lit
-from pairs.ir.loops import ParticleFor
+from pairs.ir.loops import For, ParticleFor
 from pairs.ir.operators import Operators
 from pairs.ir.operator_class import OperatorClass
+from pairs.ir.properties import ContactProperty
 from pairs.ir.scalars import ScalarOp
 from pairs.ir.types import Types
 from pairs.mapping.keywords import Keywords
@@ -228,7 +229,36 @@ class BuildParticleIR(ast.NodeVisitor):
 
     def visit_Subscript(self, node):
         #print(ast.dump(node))
-        return self.visit(node.value)[self.visit(node.slice)]
+        value = self.visit(node.value)
+        index = self.visit(node.slice)
+
+        if isinstance(value, ContactProperty) and self.sim.neighbor_lists is None:
+            i = index[0]
+            j = index[1]
+
+            if '__contact_id__' not in self.ctx_symbols:
+                particle_uid = self.sim.particle_uid
+                contact_lists = self.sim._contact_history.contact_lists
+                num_contacts = self.sim._contact_history.num_contacts
+                contact_id = self.sim.add_temp_var(-1)
+
+                for k in For(self.sim, 0, num_contacts[i]):
+                    for _ in Filter(self.sim, ScalarOp.cmp(contact_lists[i][k], particle_uid[j])):
+                        Assign(self.sim, contact_id, k)
+
+                for _ in Filter(self.sim, ScalarOp.cmp(contact_id, -1)):
+                    last_contact = num_contacts[i]
+                    Assign(self.sim, contact_id, last_contact)
+                    Assign(self.sim, num_contacts[i], last_contact + 1)
+
+                    for contact_prop in self.sim.contact_properties:
+                        Assign(self.sim, contact_prop[i, last_contact], contact_prop.default())
+
+                self.ctx_symbols.update({'__contact_id__': contact_id})
+
+            return value[i, self.ctx_symbols['__contact_id__']]
+
+        return value[index]
 
     def visit_Tuple(self, node):
         #print(ast.dump(node))
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 682d6db..85b333f 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -324,7 +324,9 @@ class Simulation:
             timestep_procedures.append((BuildNeighborLists(self, self.neighbor_lists), self.reneighbor_frequency))
 
         if self._use_contact_history:
-            timestep_procedures.append((BuildContactHistory(self, self._contact_history, self.cell_lists), self.reneighbor_frequency))
+            if self.neighbor_lists is not None:
+                timestep_procedures.append((BuildContactHistory(self, self._contact_history, self.cell_lists), self.reneighbor_frequency))
+
             timestep_procedures.append(ResetContactHistoryUsageStatus(self, self._contact_history))
 
         timestep_procedures += [
-- 
GitLab