diff --git a/data/planes.input b/data/planes.input
index 6a52e89546296d49a36239f9991f8bc22497b0ed..5de074bf8a14746a8f30ad3e3f612ff34a0c70c9 100644
--- a/data/planes.input
+++ b/data/planes.input
@@ -1,2 +1,2 @@
-0,1,0,0,0,0,0,1,13
-0,1,0.8,0.015,0.2,0,0,-1,13
+0,1,0.0,0.0,0.0,0.0,0.0,1.0,13
+0,1,0.8,0.015,0.2,0.0,0.0,-1.0,13
diff --git a/examples/dem.py b/examples/dem.py
index a8e78b882203923787cd61021c93507a73364a7f..07b68d75f0cb1d82c8dd58dd82b7908c6c1bb6aa 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -12,7 +12,7 @@ def update_mass_and_inertia(i):
inv_inertia[i] = inversed(diagonal_matrix(0.4 * mass[i] * radius[i] * radius[i]))
else:
- #mass[i] = infinity
+ mass[i] = infinity
inv_inertia[i] = 0.0
@@ -112,7 +112,7 @@ restitutionCoefficient = 0.1
collisionTime_SI = 5e-4
poissonsRatio = 0.22
timeSteps = 10000
-visSpacing = 1
+visSpacing = 100
denseBottomLayer = False
bottomLayerOffsetFactor = 1.0
kappa = 2.0 * (1.0 - poissonsRatio) / (2.0 - poissonsRatio) # from Thornton et al
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index d38b06134d61acadef76ee93b5c77f1780b707a0..eae9daa622d315547357875fdde27464c7ec824b 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -206,12 +206,12 @@ class CGen:
for array_access in kernel.array_accesses():
type_kw = Types.c_keyword(array_access.type())
- decl = f"{type_kw} a{array_access.id()}"
+ decl = f"{type_kw} {array_access.name()}"
kernel_params += f", {decl}"
for scalar_op in kernel.scalar_ops():
type_kw = Types.c_keyword(scalar_op.type())
- decl = f"{type_kw} e{scalar_op.id()}"
+ decl = f"{type_kw} {scalar_op.name()}"
kernel_params += f", {decl}"
self.print(f"__global__ void {kernel.name}({kernel_params}) {{")
@@ -268,7 +268,7 @@ class CGen:
array_name = self.generate_expression(array_access.array)
tkw = Types.c_keyword(array_access.type())
acc_index = self.generate_expression(array_access.flat_index)
- acc_ref = f"a{array_access.id()}"
+ acc_ref = array_access.name()
self.print(f"const {tkw} {acc_ref} = {array_name}[{acc_index}];")
if isinstance(ast_node.elem, AtomicAdd):
@@ -276,7 +276,7 @@ class CGen:
elem = self.generate_expression(atomic_add.elem)
value = self.generate_expression(atomic_add.value)
tkw = Types.c_keyword(atomic_add.type())
- acc_ref = f"atm_add{atomic_add.id()}"
+ acc_ref = atomic_add.name()
prefix = "" if ast_node.elem.device_flag else "host_"
if atomic_add.check_for_resize():
@@ -290,7 +290,7 @@ class CGen:
contact_prop_access = ast_node.elem
contact_prop = contact_prop_access.contact_prop
prop_name = self.generate_expression(contact_prop)
- acc_ref = f"cp{contact_prop_access.id()}"
+ acc_ref = contact_prop_access.name()
if not contact_prop_access.is_scalar():
for dim in contact_prop_access.indexes_to_generate():
@@ -306,7 +306,7 @@ class CGen:
feature_prop_access = ast_node.elem
feature_prop = feature_prop_access.feature_prop
prop_name = self.generate_expression(feature_prop)
- acc_ref = f"f{feature_prop_access.id()}"
+ acc_ref = feature_prop_access.name()
if not feature_prop_access.is_scalar():
for dim in feature_prop_access.indexes_to_generate():
@@ -321,7 +321,7 @@ class CGen:
if isinstance(ast_node.elem, PropertyAccess):
prop_access = ast_node.elem
prop_name = self.generate_expression(prop_access.prop)
- acc_ref = f"p{prop_access.id()}"
+ acc_ref = prop_access.name()
if not prop_access.is_scalar():
for dim in prop_access.indexes_to_generate():
@@ -336,7 +336,7 @@ class CGen:
quaternion = ast_node.elem
for i in quaternion.indexes_to_generate():
expr = self.generate_expression(quaternion.get_value(i))
- self.print(f"const double q{quaternion.id()}_{i} = {expr};")
+ self.print(f"const double {quaternion.name()}_{i} = {expr};")
if isinstance(ast_node.elem, QuaternionOp):
quat_op = ast_node.elem
@@ -346,9 +346,9 @@ class CGen:
operator = quat_op.operator()
if operator.is_unary():
- self.print(f"const double e{quat_op.id()}_{dim} = {operator.symbol()}({lhs});")
+ self.print(f"const double {quat_op.name()}_{dim} = {operator.symbol()}({lhs});")
else:
- self.print(f"const double e{quat_op.id()}_{dim} = {lhs} {operator.symbol()} {rhs};")
+ self.print(f"const double {quat_op.name()}_{dim} = {lhs} {operator.symbol()} {rhs};")
if isinstance(ast_node.elem, ScalarOp):
scalar_op = ast_node.elem
@@ -359,13 +359,13 @@ class CGen:
tkw = Types.c_keyword(scalar_op.type())
if operator.is_unary():
- self.print(f"const {tkw} e{scalar_op.id()} = {operator.symbol()}({lhs});")
+ self.print(f"const {tkw} {scalar_op.name()} = {operator.symbol()}({lhs});")
else:
- self.print(f"const {tkw} e{scalar_op.id()} = {lhs} {operator.symbol()} {rhs};")
+ self.print(f"const {tkw} {scalar_op.name()} = {lhs} {operator.symbol()} {rhs};")
if isinstance(ast_node.elem, Select):
select = ast_node.elem
- acc_ref = f"s{select.id()}"
+ acc_ref = select.name()
if not select.is_scalar():
for dim in select.indexes_to_generate():
@@ -382,16 +382,15 @@ class CGen:
if isinstance(ast_node.elem, MathFunction):
math_func = ast_node.elem
- acc_ref = f"mf{math_func.id()}"
params = ", ".join([str(self.generate_expression(p)) for p in math_func.parameters()])
tkw = Types.c_keyword(math_func.type())
- self.print(f"const {tkw} {acc_ref} = {math_func.function_name()}({params});")
+ self.print(f"const {tkw} {math_func.name()} = {math_func.function_name()}({params});")
if isinstance(ast_node.elem, Matrix):
matrix = ast_node.elem
for i in matrix.indexes_to_generate():
expr = self.generate_expression(matrix.get_value(i))
- self.print(f"const double m{matrix.id()}_{i} = {expr};")
+ self.print(f"const double {matrix.name()}_{i} = {expr};")
if isinstance(ast_node.elem, MatrixOp):
matrix_op = ast_node.elem
@@ -401,15 +400,15 @@ class CGen:
operator = vector_op.operator()
if operator.is_unary():
- self.print(f"const double e{matrix_op.id()}_{dim} = {operator.symbol()}({lhs});")
+ self.print(f"const double {matrix_op.name()}_{dim} = {operator.symbol()}({lhs});")
else:
- self.print(f"const double e{matrix_op.id()}_{dim} = {lhs} {operator.symbol()} {rhs};")
+ self.print(f"const double {matrix_op.name()}_{dim} = {lhs} {operator.symbol()} {rhs};")
if isinstance(ast_node.elem, Vector):
vector = ast_node.elem
for dim in vector.indexes_to_generate():
expr = self.generate_expression(vector.get_value(dim))
- self.print(f"const double v{vector.id()}_{dim} = {expr};")
+ self.print(f"const double {vector.name()}_{dim} = {expr};")
if isinstance(ast_node.elem, VectorOp):
vector_op = ast_node.elem
@@ -419,9 +418,9 @@ class CGen:
operator = vector_op.operator()
if operator.is_unary():
- self.print(f"const double e{vector_op.id()}_{dim} = {operator.symbol()}({lhs});")
+ self.print(f"const double {vector_op.name()}_{dim} = {operator.symbol()}({lhs});")
else:
- self.print(f"const double e{vector_op.id()}_{dim} = {lhs} {operator.symbol()} {rhs};")
+ self.print(f"const double {vector_op.name()}_{dim} = {lhs} {operator.symbol()} {rhs};")
if isinstance(ast_node, Branch):
cond = self.generate_expression(ast_node.cond)
@@ -745,10 +744,10 @@ class CGen:
acc_index = self.generate_expression(ast_node.flat_index)
return f"{array_name}[{acc_index}]"
- return f"a{ast_node.id()}"
+ return f"{ast_node.name()}"
if isinstance(ast_node, AtomicAdd):
- return f"atm_add{ast_node.id()}"
+ return f"{ast_node.name()}"
if isinstance(ast_node, ScalarOp):
if ast_node.inlined is True:
@@ -758,7 +757,7 @@ class CGen:
return f"({operator.symbol()}({lhs}))" if operator.is_unary() else \
f"({lhs} {operator.symbol()} {rhs})"
- return f"e{ast_node.id()}"
+ return f"{ast_node.name()}"
if isinstance(ast_node, Call):
extra_params = []
@@ -797,7 +796,7 @@ class CGen:
if isinstance(ast_node, Iter):
assert mem is False, "Iterator is not lvalue!"
- return f"i{ast_node.id()}"
+ return f"{ast_node.name()}"
if isinstance(ast_node, Lit):
assert mem is False, "Literal is not lvalue!"
@@ -820,7 +819,7 @@ class CGen:
params = ", ".join([str(self.generate_expression(p)) for p in ast_node.parameters()])
return f"{ast_node.function_name()}({params})"
- return f"mf{ast_node.id()}"
+ return f"{ast_node.name()}"
if isinstance(ast_node, Property):
return ast_node.name()
@@ -837,7 +836,7 @@ class CGen:
return f"{prop_name}[{index_expr}]"
- return f"p{ast_node.id()}" + (f"_{index}" if not ast_node.is_scalar() else "")
+ return f"{ast_node.name()}" + (f"_{index}" if not ast_node.is_scalar() else "")
if isinstance(ast_node, ContactPropertyAccess):
assert ast_node.is_scalar() or index is not None, \
@@ -851,7 +850,7 @@ class CGen:
return f"{prop_name}[{index_expr}]"
- return f"cp{ast_node.id()}" + (f"_{index}" if not ast_node.is_scalar() else "")
+ return f"{ast_node.name()}" + (f"_{index}" if not ast_node.is_scalar() else "")
if isinstance(ast_node, FeaturePropertyAccess):
assert ast_node.is_scalar() or index is not None, \
@@ -865,7 +864,7 @@ class CGen:
return f"{feature_name}[{index_expr}]"
- return f"f{ast_node.id()}" + (f"_{index}" if not ast_node.is_scalar() else "")
+ return f"{ast_node.name()}" + (f"_{index}" if not ast_node.is_scalar() else "")
if isinstance(ast_node, ParticleAttributeList):
tid = CGen.temp_id
@@ -892,9 +891,9 @@ class CGen:
if not ast_node.is_scalar():
assert index is not None, "Index must be set for non-scalar reference."
- return f"s{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
- return f"s{ast_node.id()}"
+ return f"{ast_node.name()}"
if isinstance(ast_node, Var):
return ast_node.name() if ast_node.is_scalar() else f"{ast_node.name()}_{index}"
@@ -910,27 +909,27 @@ class CGen:
if isinstance(ast_node, Vector):
assert index is not None, "Index must be set for vector."
- return f"v{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
if isinstance(ast_node, Matrix):
assert index is not None, "Index must be set for matrix."
- return f"m{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
if isinstance(ast_node, Quaternion):
assert index is not None, "Index must be set for quaternion."
- return f"q{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
if isinstance(ast_node, VectorOp):
assert index is not None, "Index must be set for vector operation."
- return f"e{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
if isinstance(ast_node, MatrixOp):
assert index is not None, "Index must be set for matrix operation."
- return f"e{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
if isinstance(ast_node, QuaternionOp):
assert index is not None, "Index must be set for quaternion operation."
- return f"e{ast_node.id()}_{index}"
+ return f"{ast_node.name()}_{index}"
if isinstance(ast_node, ZeroVector):
return "0.0"
diff --git a/src/pairs/ir/arrays.py b/src/pairs/ir/arrays.py
index 445fb96954bac8fce908327b75860b288a687538..b7e1e4245630323a983aa1bba36115162e7a783b 100644
--- a/src/pairs/ir/arrays.py
+++ b/src/pairs/ir/arrays.py
@@ -190,6 +190,9 @@ class ArrayAccess(ASTTerm):
def id(self):
return self.acc_id
+ def name(self):
+ return f"arr_acc{self.id()}" + self.label_suffix()
+
def type(self):
return self.array.type()
diff --git a/src/pairs/ir/ast_term.py b/src/pairs/ir/ast_term.py
index 778bdfc01da3d1d0476cd77e324a94e693ed6463..4c2c7b407a5357650f152b8f49dcd07c4b458a55 100644
--- a/src/pairs/ir/ast_term.py
+++ b/src/pairs/ir/ast_term.py
@@ -8,6 +8,7 @@ class ASTTerm(ASTNode):
super().__init__(sim)
self._class_type = class_type
self._indexes_to_generate = set()
+ self._label = None
def __add__(self, other):
return self._class_type(self.sim, self, other, Operators.Add)
@@ -84,6 +85,12 @@ class ASTTerm(ASTNode):
def is_quaternion(self):
return self.type() == Types.Quaternion
+ def set_label(self, label):
+ self._label = label
+
+ def label_suffix(self):
+ return "" if self._label is None else f"_{self._label}"
+
def indexes_to_generate(self):
return self._indexes_to_generate
diff --git a/src/pairs/ir/atomic.py b/src/pairs/ir/atomic.py
index 4d5c8ce9e27575b4e70292b12acb034d004e6ade..88c54b6180f3fce3d24aa0a2678c546667cf377f 100644
--- a/src/pairs/ir/atomic.py
+++ b/src/pairs/ir/atomic.py
@@ -32,6 +32,9 @@ class AtomicAdd(ASTTerm):
def id(self):
return self.atomic_add_id
+ def name(self):
+ return f"atm_add{self.id()}" + self.label_suffix()
+
def type(self):
return self.elem.type()
diff --git a/src/pairs/ir/features.py b/src/pairs/ir/features.py
index 5a2a29270a5ce788b49b1c7f2ae9a40853bebbb1..7067073b67fddfbc034a2230c621e0c4f65f1fee 100644
--- a/src/pairs/ir/features.py
+++ b/src/pairs/ir/features.py
@@ -185,6 +185,9 @@ class FeaturePropertyAccess(ASTTerm):
def id(self):
return self.acc_id
+ def name(self):
+ return f"feat_prop_acc{self.id()}" + self.label_suffix()
+
def type(self):
return self.feature_prop.type()
diff --git a/src/pairs/ir/math.py b/src/pairs/ir/math.py
index b1cab83a481d578e719445b2c4da64fa9a6141fa..9dde4d2bedcbd7e7951e77d1bad5f5cdf3469101 100644
--- a/src/pairs/ir/math.py
+++ b/src/pairs/ir/math.py
@@ -23,6 +23,9 @@ class MathFunction(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"mf{self.id()}" + self.label_suffix()
+
def function_name(self):
return "undefined"
diff --git a/src/pairs/ir/matrices.py b/src/pairs/ir/matrices.py
index 88baf50dd034287e3791682e8521c4cbaf74d41b..aff53add5cc3db80ef66b098bed1b4e622cf66b7 100644
--- a/src/pairs/ir/matrices.py
+++ b/src/pairs/ir/matrices.py
@@ -35,6 +35,9 @@ class MatrixOp(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"mat_op{self.id()}" + self.label_suffix()
+
def type(self):
return Types.Matrix
@@ -100,6 +103,9 @@ class Matrix(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"mat{self.id()}" + self.label_suffix()
+
def type(self):
return Types.Matrix
diff --git a/src/pairs/ir/properties.py b/src/pairs/ir/properties.py
index 396510df2d850afdb535ce517003785b1040947c..b4867f8be15acd56e5468f0c0f830163e1a5c4cd 100644
--- a/src/pairs/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -30,6 +30,9 @@ class Properties:
def volatiles(self):
return [p for p in self.props if p.volatile is True]
+ def non_volatiles(self):
+ return [p for p in self.props if p.volatile is False]
+
def nprops(self):
return len(self.props)
@@ -135,6 +138,9 @@ class PropertyAccess(ASTTerm):
def id(self):
return self.acc_id
+ def name(self):
+ return f"prop_acc{self.id()}" + self.label_suffix()
+
def type(self):
return self.prop.type()
@@ -298,6 +304,9 @@ class ContactPropertyAccess(ASTTerm):
def id(self):
return self.acc_id
+ def name(self):
+ return f"cont_prop_acc{self.id()}" + self.label_suffix()
+
def type(self):
return self.contact_prop.type()
diff --git a/src/pairs/ir/quaternions.py b/src/pairs/ir/quaternions.py
index ef9a7e02f27a375982a87ddf0eb2b24bbb77c6cf..5406973f160ab306a97fd52fc443dc140ba0a383 100644
--- a/src/pairs/ir/quaternions.py
+++ b/src/pairs/ir/quaternions.py
@@ -35,6 +35,9 @@ class QuaternionOp(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"quat_op{self.id()}" + self.label_suffix()
+
def type(self):
return Types.Quaternion
@@ -100,6 +103,9 @@ class Quaternion(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"quat{self.id()}" + self.label_suffix()
+
def type(self):
return Types.Quaternion
diff --git a/src/pairs/ir/scalars.py b/src/pairs/ir/scalars.py
index 1e33274122c07f25dde5114e9f10f11dda440c85..12c2462672ab50b3ade5f7b20fc8db36714666eb 100644
--- a/src/pairs/ir/scalars.py
+++ b/src/pairs/ir/scalars.py
@@ -116,6 +116,9 @@ class ScalarOp(ASTTerm):
def id(self):
return self.scalar_op_id
+ def name(self):
+ return f"sca_op{self.id()}" + self.label_suffix()
+
def type(self):
return self.scalar_op_type
diff --git a/src/pairs/ir/select.py b/src/pairs/ir/select.py
index 9bc43e10d849c0aa60ea2c288a4adcd9f9e01c81..e920f00b460b3145335556d412af7afb2a65f586 100644
--- a/src/pairs/ir/select.py
+++ b/src/pairs/ir/select.py
@@ -36,6 +36,9 @@ class Select(ASTTerm):
def id(self):
return self.select_id
+ def name(self):
+ return f"sel{self.id()}" + self.label_suffix()
+
def type(self):
return self.expr_if.type()
diff --git a/src/pairs/ir/vectors.py b/src/pairs/ir/vectors.py
index 825cb9716f91718794d33a2352db4dc2266cd97b..53e5f57fd6dbe7a6fce847aa5ecfa4bd37938afd 100644
--- a/src/pairs/ir/vectors.py
+++ b/src/pairs/ir/vectors.py
@@ -35,6 +35,9 @@ class VectorOp(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"vec_op{self.id()}" + self.label_suffix()
+
def type(self):
return Types.Vector
@@ -108,6 +111,9 @@ class Vector(ASTTerm):
def id(self):
return self._id
+ def name(self):
+ return f"vec{self.id()}" + self.label_suffix()
+
def type(self):
return Types.Vector
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index 62c7c2bfb535fb46a460fc95a7f6fc97bd68a227..41152b8fcb5a3bc0044d107423d17b9e311096cb 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -95,6 +95,7 @@ class BuildParticleIR(ast.NodeVisitor):
if isinstance(lhs, UndefinedSymbol):
self.add_symbols({lhs.symbol_id: rhs})
+ rhs.set_label(lhs.symbol_id)
else:
Assign(self.sim, lhs, rhs)
@@ -106,6 +107,7 @@ class BuildParticleIR(ast.NodeVisitor):
if isinstance(lhs, UndefinedSymbol):
self.add_symbols({lhs.symbol_id: bin_op})
+ rhs.set_label(lhs.symbol_id)
else:
Assign(self.sim, lhs, bin_op)
diff --git a/src/pairs/sim/comm.py b/src/pairs/sim/comm.py
index d945a404d335ea0b0294b025eb74258290328845..06186060b6bdeaa30aaad60511d1b2f7fc2df505 100644
--- a/src/pairs/sim/comm.py
+++ b/src/pairs/sim/comm.py
@@ -46,9 +46,10 @@ class Comm:
@pairs_inline
def borders(self):
- # Every property that is constant across timesteps and have neighbor accesses during any
- # interaction kernel (i.e. property[j] in force calculation kernel)
- prop_list = [self.sim.property(p) for p in ['mass', 'position', 'linear_velocity', 'angular_velocity', 'flags']]
+ # Every property that has neighbor accesses during any interaction kernel (i.e. property[j]
+ # exists in any force calculation kernel)
+ # We ignore normal because there should be no halfspace ghosts
+ prop_list = [self.sim.property(p) for p in ['mass', 'radius', 'position', 'linear_velocity', 'angular_velocity', 'shape']]
Assign(self.sim, self.nsend_all, 0)
Assign(self.sim, self.sim.nghost, 0)
@@ -64,7 +65,7 @@ class Comm:
@pairs_inline
def exchange(self):
# Every property except volatiles
- prop_list = [self.sim.property(p) for p in ['mass', 'position', 'linear_velocity', 'shape', 'flags']]
+ prop_list = self.sim.properties.non_volatiles()
for step in range(self.dom_part.number_of_steps()):
Assign(self.sim, self.nsend_all, 0)
Assign(self.sim, self.sim.nghost, 0)