diff --git a/examples/lj.py b/examples/lj.py
index 98b2b58c2b18acc130afbeb4f04b84951c9160ac..8f2fe40fe361b9b25da4cbeab176340158c5d846 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -3,9 +3,9 @@ import sys
def lj(i, j):
- sr2 = 1.0 / rsq
+ sr2 = 1.0 / squared_distance(i, j)
sr6 = sr2 * sr2 * sr2 * sigma6[i, j]
- force[i] += delta * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j]
+ force[i] += delta(i, j) * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j]
def euler(i):
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index dbe69bbcac452ca698104ce3b2140e0abbfc582a..0accd1c9a56ae4fd644eb171a5c4503d93e28a9f 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -17,7 +17,7 @@ from pairs.ir.math import Ceil, Sqrt
from pairs.ir.memory import Malloc, Realloc
from pairs.ir.module import ModuleCall
from pairs.ir.particle_attributes import ParticleAttributeList
-from pairs.ir.properties import Property, PropertyAccess, RegisterProperty, ReallocProperty
+from pairs.ir.properties import Property, PropertyAccess, RegisterProperty, ReallocProperty, ContactPropertyAccess, RegisterContactProperty
from pairs.ir.select import Select
from pairs.ir.sizeof import Sizeof
from pairs.ir.types import Types
diff --git a/src/pairs/ir/atomic.py b/src/pairs/ir/atomic.py
index 9de1c9eb43def2ad4e21be2352a6bf4578654a7f..2bd54ae00f0430608c8462511166e440d210ed7d 100644
--- a/src/pairs/ir/atomic.py
+++ b/src/pairs/ir/atomic.py
@@ -19,7 +19,7 @@ class AtomicAdd(ASTTerm):
self.device_flag = False
def __str__(self):
- return f"AtomicAdd<{self.elem, self.val}>"
+ return f"AtomicAdd<{self.elem, self.value}>"
def add_resize_check(self, resize, capacity):
self.resize = BinOp.inline(resize)
diff --git a/src/pairs/ir/bin_op.py b/src/pairs/ir/bin_op.py
index 66673b285f0c45588dd92e2c63bcb338e3134eda..4b3b32b73a778823d58871cae654772d0e81345d 100644
--- a/src/pairs/ir/bin_op.py
+++ b/src/pairs/ir/bin_op.py
@@ -52,7 +52,7 @@ class BinOp(VectorExpression):
def __str__(self):
a = self.lhs.id() if isinstance(self.lhs, BinOp) else self.lhs
b = self.rhs.id() if isinstance(self.rhs, BinOp) else self.rhs
- return f"BinOp<{a} {self.op} {b}>"
+ return f"BinOp<{a} {self.op.symbol()} {b}>"
def copy(self):
return BinOp(self.sim, self.lhs.copy(), self.rhs.copy(), self.op, self.mem)
diff --git a/src/pairs/ir/properties.py b/src/pairs/ir/properties.py
index 35a79d5de704e3cd5cee75dcaa91fbd15ac5bd44..9bcf2ec93eff17bdd734657e93a8c43399cea641 100644
--- a/src/pairs/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -11,7 +11,6 @@ class Properties:
def __init__(self, sim):
self.sim = sim
self.props = []
- self.capacities = []
self.defs = {}
def add(self, p_name, p_type, p_value, p_volatile, p_layout=Layouts.AoS):
@@ -20,12 +19,6 @@ class Properties:
self.defs[p_name] = p_value
return p
- def add_capacity(self, var):
- self.capacities.append(var)
-
- def is_capacity(self, var):
- return var in self.capacities
-
def defaults(self):
return self.defs
@@ -207,6 +200,9 @@ class ContactProperties:
return None
+ def empty(self):
+ return len(self.contact_properties) == 0
+
def __iter__(self):
yield from self.contact_properties
diff --git a/src/pairs/sim/contact_history.py b/src/pairs/sim/contact_history.py
index 63727028499e71cc9ed3f3380ab995bec6a5e5ec..b626a6c3c85fa4d3ce7423a0c9b0b926be250eec 100644
--- a/src/pairs/sim/contact_history.py
+++ b/src/pairs/sim/contact_history.py
@@ -1,3 +1,4 @@
+from pairs.ir.bin_op import BinOp
from pairs.ir.block import pairs_device_block
from pairs.ir.branches import Branch, Filter
from pairs.ir.loops import ParticleFor, For
@@ -37,16 +38,23 @@ class BuildContactHistory(Lowerable):
j = neigh.particle_index()
neighbor_contact.set(-1)
for k in For(self.sim, 0, num_contacts[i]):
- if contact_lists[i][k] == j:
+ for _ in Filter(self.sim, BinOp.cmp(contact_lists[i][k], j)):
neighbor_contact.set(k)
for _ in Filter(self.sim, neighbor_contact >= 0):
contact_lists[i][k].set(contact_lists[i][last_contact_id])
for contact_prop in self.sim.contact_properties:
- tmp = self.sim.add_temp_var(contact_prop[i, last_contact_id])
- contact_prop[i, last_contact_id].set(contact_prop[i, k])
- contact_prop[i, k].set(tmp)
+ if contact_prop.type() == Types.Vector:
+ for d in range(0, self.sim.ndims()):
+ tmp = self.sim.add_temp_var(contact_prop[i, last_contact_id][d])
+ contact_prop[i, last_contact_id][d].set(contact_prop[i, k][d])
+ contact_prop[i, k].set(tmp)
+
+ else:
+ tmp = self.sim.add_temp_var(contact_prop[i, last_contact_id])
+ contact_prop[i, last_contact_id].set(contact_prop[i, k])
+ contact_prop[i, k].set(tmp)
contact_lists[i][last_contact_id].set(j)
@@ -57,12 +65,16 @@ class BuildContactHistory(Lowerable):
j = neigh.particle_index()
neighbor_contact.set(-1)
for k in For(self.sim, 0, num_contacts[i]):
- if contact_lists[i][k] == j:
+ for _ in Filter(self.sim, BinOp.cmp(contact_lists[i][k], j)):
neighbor_contact.set(k)
for _ in Filter(self.sim, neighbor_contact < 0):
- for contact_prop in self.sim.contact_properties():
- contact_prop[i, last_contact_id].set(contact_prop.default())
+ for contact_prop in self.sim.contact_properties:
+ if contact_prop.type() == Types.Vector:
+ for d in range(0, self.sim.ndims()):
+ contact_prop[i, last_contact_id][d].set(contact_prop.default()[d])
+ else:
+ contact_prop[i, last_contact_id].set(contact_prop.default())
contact_lists[i][last_contact_id].set(j)
diff --git a/src/pairs/sim/properties.py b/src/pairs/sim/properties.py
index 37b28d26854b356dbe0056c40a9711f0e9613009..ec4d7dd28acf85f210239e848326454dc836ba33 100644
--- a/src/pairs/sim/properties.py
+++ b/src/pairs/sim/properties.py
@@ -16,13 +16,12 @@ class PropertiesAlloc(FinalLowerable):
@pairs_inline
def lower(self):
- capacity = sum(self.sim.properties.capacities)
for p in self.sim.properties.all():
sizes = []
if Types.is_real(p.type()) or Types.is_integer(p.type()):
- sizes = [capacity]
+ sizes = [self.sim.particle_capacity]
elif p.type() == Types.Vector:
- sizes = [capacity, self.sim.ndims()]
+ sizes = [self.sim.particle_capacity, self.sim.ndims()]
else:
raise Exception("Invalid property type!")
@@ -38,13 +37,12 @@ class ContactPropertiesAlloc(FinalLowerable):
@pairs_inline
def lower(self):
- capacity = sum(self.sim.contact_properties.capacities)
- for p in self.sim.contact_properties.all():
+ for p in self.sim.contact_properties:
sizes = []
if Types.is_real(p.type()) or Types.is_integer(p.type()):
- sizes = [capacity]
+ sizes = [self.sim.particle_capacity * self.sim.neighbor_capacity]
elif p.type() == Types.Vector:
- sizes = [capacity, self.sim.ndims()]
+ sizes = [self.sim.particle_capacity * self.sim.neighbor_capacity, self.sim.ndims()]
else:
raise Exception("Invalid contact property type!")
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 1a3cbc1f7ec1754f7915d7255256934744d24a6e..4b44eef928cc5893aab6512e151d9078319a947f 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -31,7 +31,7 @@ from pairs.transformations import Transformations
class Simulation:
- def __init__(self, code_gen, dims=3, timesteps=100, particle_capacity=1000000):
+ def __init__(self, code_gen, dims=3, timesteps=100):
self.code_gen = code_gen
self.code_gen.assign_simulation(self)
self.position_prop = None
@@ -41,8 +41,8 @@ class Simulation:
self.features = Features(self)
self.feature_properties = FeatureProperties(self)
self.contact_properties = ContactProperties(self)
- self.particle_capacity = self.add_var('particle_capacity', Types.Int32, particle_capacity)
- self.neighbor_capacity = self.add_var('neighbor_capacity', Types.Int32, particle_capacity)
+ self.particle_capacity = self.add_var('particle_capacity', Types.Int32, 200000)
+ self.neighbor_capacity = self.add_var('neighbor_capacity', Types.Int32, 100)
self.nlocal = self.add_var('nlocal', Types.Int32)
self.nghost = self.add_var('nghost', Types.Int32)
self.resizes = self.add_array('resizes', 3, Types.Int32, arr_sync=False)
@@ -71,7 +71,6 @@ class Simulation:
self._target = None
self._dom_part = DimensionRanges(self)
self.nparticles = self.nlocal + self.nghost
- self.properties.add_capacity(self.particle_capacity)
def add_module(self, module):
assert isinstance(module, Module), "add_module(): Given parameter is not of type Module!"
@@ -279,18 +278,24 @@ class Simulation:
def generate(self):
assert self._target is not None, "Target not specified!"
comm = Comm(self, self._dom_part)
- contact_history = ContactHistory(self.cell_lists)
- timestep = Timestep(self, self.ntimesteps, [
+ timestep_procedures = [
(comm.exchange(), 20),
(comm.borders(), comm.synchronize(), 20),
(CellListsBuild(self, self.cell_lists), 20),
(NeighborListsBuild(self, self.neighbor_lists), 20),
- (BuildContactHistory(self, contact_history), 20),
+ ]
+
+ if not self.contact_properties.empty():
+ contact_history = ContactHistory(self.cell_lists)
+ timestep_procedures.append((BuildContactHistory(self, contact_history), 20))
+
+ timestep_procedures += [
PropertiesResetVolatile(self),
self.functions
- ])
+ ]
+ timestep = Timestep(self, self.ntimesteps, timestep_procedures)
self.enter(timestep.block)
timestep.add(VTKWrite(self, self.vtk_file, timestep.timestep() + 1))
self.leave()
diff --git a/src/pairs/transformations/modules.py b/src/pairs/transformations/modules.py
index 387b8f6b16ed7db0b294d89ee74dba94178f6abf..294f46c7dc5335f05e244b3e5cd795b90ae39955 100644
--- a/src/pairs/transformations/modules.py
+++ b/src/pairs/transformations/modules.py
@@ -166,9 +166,9 @@ class ReplaceModulesByCalls(Mutator):
branch_cond = cond if branch_cond is None else BinOp.or_op(cond, branch_cond)
props_realloc = []
- if properties.is_capacity(capacity):
+ if capacity == sim.particle_capacity:
for p in properties.all():
- new_capacity = sum(properties.capacities)
+ new_capacity = sim.particle_capacity
sizes = [new_capacity, sim.ndims()] if p.type() == Types.Vector else [new_capacity]
props_realloc += [ReallocProperty(sim, p, sizes)]