From 6b2c436e83470d452fc86bdbe0eea033bc66145c Mon Sep 17 00:00:00 2001
From: Behzad Safaei <iwia103h@alex2.nhr.fau.de>
Date: Fri, 7 Feb 2025 03:23:00 +0100
Subject: [PATCH] Add feature props to accessor, fix issue with non scalar
types
---
examples/modular/sd_sample_3_GPU.cu | 20 ++++-
examples/modular/spring_dashpot.py | 18 ++--
src/pairs/code_gen/accessor.py | 120 +++++++++++++++++++++++++++
src/pairs/code_gen/cgen.py | 9 +-
src/pairs/ir/device.py | 22 +++++
src/pairs/ir/features.py | 6 +-
src/pairs/sim/simulation.py | 10 ++-
src/pairs/transformations/devices.py | 5 +-
8 files changed, 190 insertions(+), 20 deletions(-)
diff --git a/examples/modular/sd_sample_3_GPU.cu b/examples/modular/sd_sample_3_GPU.cu
index 964ede0..6861c08 100644
--- a/examples/modular/sd_sample_3_GPU.cu
+++ b/examples/modular/sd_sample_3_GPU.cu
@@ -24,6 +24,20 @@ __global__ void change_gravitational_force(PairsAccessor ac, int idx){
printf("Force [from device] after setting = (%f, %f, %f) \n", ac.getForce(idx)[0], ac.getForce(idx)[1], ac.getForce(idx)[2]);
}
+void set_feature_properties(std::shared_ptr<PairsAccessor> &ac){
+ ac->setTypeStiffness(0,0, 0);
+ ac->setTypeStiffness(0,1, 1000);
+ ac->setTypeStiffness(1,0, 1000);
+ ac->setTypeStiffness(1,1, 3000);
+ ac->syncTypeStiffness();
+
+ ac->setTypeDampingNorm(0,0, 0);
+ ac->setTypeDampingNorm(0,1, 20);
+ ac->setTypeDampingNorm(1,0, 20);
+ ac->setTypeDampingNorm(1,1, 10);
+ ac->syncTypeDampingNorm();
+}
+
int main(int argc, char **argv) {
auto pairs_sim = std::make_shared<PairsSimulation>();
@@ -41,8 +55,10 @@ int main(int argc, char **argv) {
pairs_sim->create_halfspace(1,1,1, 0, -1, 0, 0, 13);
pairs_sim->create_halfspace(1,1,1, 0, 0, -1, 0, 13);
- pairs::id_t pUid = pairs_sim->create_sphere(0.6, 0.6, 0.7, 0, 0, 0, 1000, 0.05, 0, 0);
- pairs_sim->create_sphere(0.4, 0.4, 0.76, 2, 2, 0, 1000, 0.05, 0, 0);
+ pairs::id_t pUid = pairs_sim->create_sphere(0.6, 0.6, 0.7, 0, 0, 0, 1000, 0.05, 1, 0);
+ pairs_sim->create_sphere(0.4, 0.4, 0.76, 2, 2, 0, 1000, 0.05, 1, 0);
+
+ set_feature_properties(ac);
MPI_Allreduce(MPI_IN_PLACE, &pUid, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_COMM_WORLD);
diff --git a/examples/modular/spring_dashpot.py b/examples/modular/spring_dashpot.py
index 94e5ef8..3109fa5 100644
--- a/examples/modular/spring_dashpot.py
+++ b/examples/modular/spring_dashpot.py
@@ -48,11 +48,6 @@ def euler(i):
def gravity(i):
force[i][2] = force[i][2] - mass[i] * gravity_SI
-gravity_SI = 9.81
-diameter = 100 # required for linkedCellWidth. TODO: set linkedCellWidth at runtime
-linkedCellWidth = 1.01 * diameter
-
-ntypes = 1
file_name = os.path.basename(__file__)
file_name_without_extension = os.path.splitext(file_name)[0]
@@ -76,6 +71,10 @@ elif target == 'cpu':
else:
print(f"Invalid target, use {sys.argv[0]} <cpu/gpu>")
+gravity_SI = 9.81
+diameter = 100 # required for linkedCellWidth. TODO: set linkedCellWidth at runtime
+linkedCellWidth = 1.01 * diameter
+ntypes = 2
psim.add_position('position')
psim.add_property('mass', pairs.real())
@@ -83,20 +82,17 @@ psim.add_property('linear_velocity', pairs.vector())
psim.add_property('angular_velocity', pairs.vector())
psim.add_property('force', pairs.vector(), volatile=True)
psim.add_property('torque', pairs.vector(), volatile=True)
-psim.add_property('hydrodynamic_force', pairs.vector(), reduce=True)
-psim.add_property('hydrodynamic_torque', pairs.vector(), reduce=True)
-psim.add_property('old_hydrodynamic_force', pairs.vector())
-psim.add_property('old_hydrodynamic_torque', pairs.vector())
psim.add_property('radius', pairs.real())
psim.add_property('normal', pairs.vector())
psim.add_property('inv_inertia', pairs.matrix())
psim.add_property('rotation_matrix', pairs.matrix())
psim.add_property('rotation', pairs.quaternion())
+
psim.add_feature('type', ntypes)
psim.add_feature_property('type', 'stiffness', pairs.real(), [3000 for i in range(ntypes * ntypes)])
psim.add_feature_property('type', 'damping_norm', pairs.real(), [10.0 for i in range(ntypes * ntypes)])
-psim.add_feature_property('type', 'damping_tan', pairs.real(), [0 for i in range(ntypes * ntypes)])
-psim.add_feature_property('type', 'friction', pairs.real(), [0 for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'damping_tan', pairs.real())
+psim.add_feature_property('type', 'friction', pairs.real())
psim.set_domain_partitioner(pairs.block_forest())
psim.pbc([False, False, False])
diff --git a/src/pairs/code_gen/accessor.py b/src/pairs/code_gen/accessor.py
index cd5a986..5948b15 100644
--- a/src/pairs/code_gen/accessor.py
+++ b/src/pairs/code_gen/accessor.py
@@ -38,12 +38,122 @@ class PairsAcessor:
self.set_property(p)
self.sync_property(p)
+ for fp in self.sim.feature_properties:
+ self.get_feature_property(fp)
+ self.set_feature_property(fp)
+ self.sync_feature_property(fp)
+
self.utility_funcs()
self.print.add_indent(-4)
self.print("};")
self.print("")
+ def get_fp_body(self, fp, device=False):
+ self.print.add_indent(4)
+ fp_name = fp.name()
+ f_name = fp.feature().name()
+
+ tkw = Types.c_accessor_keyword(self.sim, fp.type())
+
+ if self.target.is_gpu() and device:
+ v = f"{fp_name}_d"
+ else:
+ v = f"ps->pobj->{fp_name}"
+
+ idx = f"{fp.feature().nkinds()}*{f_name}1 + {f_name}2"
+ if Types.is_scalar(fp.type()):
+ self.print(f"return {v}[{idx}];")
+ else:
+ nelems = Types.number_of_elements(self.sim, fp.type())
+ return_values = [f"{v}[({idx})*{nelems} + {n}]" for n in range(nelems)]
+ self.print(f"return {tkw}(" + ", ".join(rv for rv in return_values) + ");")
+ self.print.add_indent(-4)
+
+
+ def get_feature_property(self, fp):
+ fp_name = fp.name()
+ tkw = Types.c_accessor_keyword(self.sim, fp.type())
+ f_name = fp.feature().name()
+
+ splitname = f_name.split('_') + fp_name.split('_')
+ funcname = ''.join(word.capitalize() for word in splitname)
+ self.print(f"{self.host_device_attr}{tkw} get{funcname}(const size_t {f_name}1, const size_t {f_name}2) const{{")
+
+ if self.target.is_gpu():
+ self.ifdef_else("__CUDA_ARCH__", self.get_fp_body, [fp, True], self.get_fp_body, [fp, False])
+ else:
+ self.get_fp_body(fp, False)
+
+ self.print("}")
+ self.print("")
+
+ def set_fp_body(self, fp, device=False):
+ self.print.add_indent(4)
+ fp_name = fp.name()
+ f_name = fp.feature().name()
+ tkw = Types.c_accessor_keyword(self.sim, fp.type())
+
+ if self.target.is_gpu() and device:
+ v = f"{fp_name}_d"
+ else:
+ v = f"ps->pobj->{fp_name}"
+
+ idx = f"{fp.feature().nkinds()}*{f_name}1 + {f_name}2"
+
+ if Types.is_scalar(fp.type()):
+ self.print(f"{v}[{idx}] = value;")
+ else:
+ nelems = Types.number_of_elements(self.sim, fp.type())
+ set_values = [f"{v}[({idx})*{nelems} + {n}] = value[{n}];" for n in range(nelems)]
+ for sv in set_values:
+ self.print(sv)
+
+ if self.target.is_gpu():
+ flag = f"*{fp_name}_device_flag_d" if device else f"{fp_name}_host_flag"
+ self.print(f"{flag} = true;")
+
+ self.print.add_indent(-4)
+
+
+ def set_feature_property(self, fp):
+ fp_name = fp.name()
+ tkw = Types.c_accessor_keyword(self.sim, fp.type())
+ f_name = fp.feature().name()
+
+ splitname = f_name.split('_') + fp_name.split('_')
+ funcname = ''.join(word.capitalize() for word in splitname)
+
+ # Feature properties can only be set from host
+ self.print(f"{self.host_attr}void set{funcname}(const size_t {f_name}1, const size_t {f_name}2, const {tkw} &value){{")
+ self.set_fp_body(fp, False)
+
+ self.print("}")
+ self.print("")
+
+ def sync_feature_property(self, fp):
+ fp_id = fp.id()
+ fp_name = fp.name()
+ f_name = fp.feature().name()
+ splitname = f_name.split('_') + fp_name.split('_')
+ funcname = ''.join(word.capitalize() for word in splitname)
+
+ self.print(f"{self.host_attr}void sync{funcname}(SyncMode sync_mode = HostAndDevice){{")
+
+ if self.target.is_gpu():
+ self.print.add_indent(4)
+ self.print(f"if ({fp_name}_host_flag) {{")
+ self.print(f" ps->pairs_runtime->copyFeaturePropertyToDevice({fp_id});")
+ self.print("}")
+ self.print("")
+
+ self.print(f"{fp_name}_host_flag = false;")
+ self.print.add_indent(-4)
+
+ self.print("}")
+ self.print("")
+
+
def member_variables(self):
self.print("PairsSimulation *ps;")
@@ -66,6 +176,16 @@ class PairsAcessor:
self.print(f"{tkw} *{pname}_device_flag_d;")
self.print(f"{tkw} {pname}_device_flag_h = false;")
self.print(f"{tkw} {pname}_host_flag = false;")
+
+ self.print("")
+ self.print("//Feature properties are global")
+
+ self.print("")
+ self.print("//Feature property flags")
+ for fp in self.sim.feature_properties:
+ fpname = fp.name()
+ tkw = Types.c_keyword(self.sim, Types.Boolean)
+ self.print(f"{tkw} {fpname}_host_flag = false;")
self.print("")
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 29e8b2f..9cfce8b 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -9,7 +9,7 @@ from pairs.ir.cast import Cast
from pairs.ir.contexts import Contexts
from pairs.ir.declaration import Decl
from pairs.ir.scalars import ScalarOp
-from pairs.ir.device import CopyArray, CopyContactProperty, CopyProperty, CopyVar, DeviceStaticRef, HostRef
+from pairs.ir.device import CopyArray, CopyContactProperty, CopyProperty, CopyFeatureProperty, CopyVar, DeviceStaticRef, HostRef
from pairs.ir.features import FeatureProperty, FeaturePropertyAccess, RegisterFeatureProperty
from pairs.ir.functions import Call
from pairs.ir.kernel import KernelLaunch
@@ -903,6 +903,13 @@ class CGen:
size = self.generate_expression(ast_node.prop().copy_size())
self.print(f"pairs_runtime->copyPropertyTo{ctx_suffix}({prop_id}, {action}, {size}); // {prop_name}")
+ if isinstance(ast_node, CopyFeatureProperty):
+ prop_id = ast_node.prop().id()
+ prop_name = ast_node.prop().name()
+ if ast_node.context() == Contexts.Device:
+ assert ast_node.action()==Actions.ReadOnly, "Feature properties can only be read from device."
+ self.print(f"pairs_runtime->copyFeaturePropertyToDevice({prop_id}); // {prop_name}")
+
if isinstance(ast_node, CopyVar):
var_name = ast_node.variable().name()
ctx_suffix = "Device" if ast_node.context() == Contexts.Device else "Host"
diff --git a/src/pairs/ir/device.py b/src/pairs/ir/device.py
index 4eddecf..0dfd6b9 100644
--- a/src/pairs/ir/device.py
+++ b/src/pairs/ir/device.py
@@ -77,7 +77,29 @@ class CopyProperty(ASTNode):
def children(self):
return [self._prop, self.sim.nghost, self.sim.nlocal]
+class CopyFeatureProperty(ASTNode):
+ def __init__(self, sim, prop, ctx, action):
+ super().__init__(sim)
+ self._prop = prop
+ self._ctx = ctx
+ self._action = action
+ self.sim.add_statement(self)
+ def __str__(self):
+ return f"CopyFeatureProperty<{self._prop}>"
+
+ def prop(self):
+ return self._prop
+
+ def context(self):
+ return self._ctx
+
+ def action(self):
+ return self._action
+
+ def children(self):
+ return [self._prop]
+
class CopyContactProperty(ASTNode):
def __init__(self, sim, prop, ctx, action):
super().__init__(sim)
diff --git a/src/pairs/ir/features.py b/src/pairs/ir/features.py
index d709336..24bb365 100644
--- a/src/pairs/ir/features.py
+++ b/src/pairs/ir/features.py
@@ -131,9 +131,7 @@ class FeatureProperty(ASTNode):
self.feature_prop_feature.nkinds()]
def array_size(self):
- nelems = self.feature_prop_feature.nkinds() * \
- Types.number_of_elements(self.sim, self.feature_prop_type)
- return nelems * nelems
+ return self.feature_prop_feature.nkinds()**2 * Types.number_of_elements(self.sim, self.feature_prop_type)
def __getitem__(self, expr):
return FeaturePropertyAccess(self.sim, self, expr)
@@ -161,7 +159,7 @@ class FeaturePropertyAccess(ASTTerm):
sizes = feature_prop.sizes()
layout = feature_prop.layout()
- for elem in range(Types.number_of_elements(feature_prop.type())):
+ for elem in range(Types.number_of_elements(self.sim, feature_prop.type())):
if layout == Layouts.AoS:
self.vector_indexes[elem] = self.index * sizes[0] + elem
elif layout == Layouts.SoA:
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 208f1ad..66d5f6d 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -217,10 +217,18 @@ class Simulation:
assert self.feature(feature_name) is None, f"Feature already defined: {feature_name}"
return self.features.add(feature_name, nkinds)
- def add_feature_property(self, feature_name, prop_name, prop_type, prop_data):
+ def add_feature_property(self, feature_name, prop_name, prop_type, prop_data=None):
feature = self.feature(feature_name)
assert feature is not None, f"Feature not found: {feature_name}"
assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
+
+ array_size = feature.nkinds()**2 * Types.number_of_elements(self, prop_type)
+
+ if not prop_data:
+ prop_data = [0 for i in range(array_size)]
+ else:
+ assert len(prop_data) == array_size, f"Incorrect array size for {prop_name}: Expected array size = {array_size}"
+
return self.feature_properties.add(feature, prop_name, prop_type, prop_data)
def add_contact_property(self, prop_name, prop_type, prop_default, layout=Layouts.AoS):
diff --git a/src/pairs/transformations/devices.py b/src/pairs/transformations/devices.py
index 21898c1..d33f30e 100644
--- a/src/pairs/transformations/devices.py
+++ b/src/pairs/transformations/devices.py
@@ -4,7 +4,7 @@ from pairs.ir.block import Block
from pairs.ir.branches import Branch, Filter
from pairs.ir.cast import Cast
from pairs.ir.contexts import Contexts
-from pairs.ir.device import CopyArray, CopyContactProperty, CopyProperty, CopyVar, DeviceStaticRef, HostRef
+from pairs.ir.device import CopyArray, CopyContactProperty, CopyProperty, CopyFeatureProperty, CopyVar, DeviceStaticRef, HostRef
from pairs.ir.functions import Call_Void
from pairs.ir.kernel import Kernel, KernelLaunch
from pairs.ir.lit import Lit
@@ -46,6 +46,9 @@ class AddDeviceCopies(Mutator):
for prop, action in s.module.properties().items():
new_stmts += [CopyProperty(s.sim, prop, copy_context, action)]
+ for fp, action in s.module.feature_properties().items():
+ new_stmts += [CopyFeatureProperty(s.sim, fp, copy_context, action)]
+
for contact_prop, action in s.module.contact_properties().items():
new_stmts += [CopyContactProperty(s.sim, contact_prop, copy_context, action)]
--
GitLab