Skip to content
Snippets Groups Projects
Commit e7f04d1a authored by Rafael Ravedutti's avatar Rafael Ravedutti
Browse files

Add option to split kernel with different types of interactions 


Signed-off-by: default avatarRafael Ravedutti <rafaelravedutti@gmail.com>
parent b2df7aa2
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/pairs/sim/interaction.py
+ 109
− 117
View file @ e7f04d1a
...@@ -5,6 +5,7 @@ from pairs.ir.block import Block, pairs_inline ...@@ -5,6 +5,7 @@ from pairs.ir.block import Block, pairs_inline
from pairs.ir.branches import Filter from pairs.ir.branches import Filter
from pairs.ir.loops import For, ParticleFor from pairs.ir.loops import For, ParticleFor
from pairs.ir.math import Sqrt from pairs.ir.math import Sqrt
from pairs.ir.select import Select
from pairs.ir.types import Types from pairs.ir.types import Types
from pairs.ir.vectors import Vector from pairs.ir.vectors import Vector
from pairs.sim.flags import Flags from pairs.sim.flags import Flags
...@@ -40,14 +41,12 @@ class Neighbor(ASTTerm): ...@@ -40,14 +41,12 @@ class Neighbor(ASTTerm):
class NeighborFor: class NeighborFor:
def number_of_cases(sim, use_cell_lists): def __init__(self, sim, particle, cell_lists, neighbor_lists=None, shapes=None):
return 2 if sim.neighbor_lists is None or use_cell_lists else 1
def __init__(self, sim, particle, cell_lists, neighbor_lists=None):
self.sim = sim self.sim = sim
self.particle = particle self.particle = particle
self.cell_lists = cell_lists self.cell_lists = cell_lists
self.neighbor_lists = neighbor_lists self.neighbor_lists = neighbor_lists
self.shapes = range(sim.max_shapes()) if shapes is None else shapes
def __str__(self): def __str__(self):
return f"NeighborFor<{self.particle}>" return f"NeighborFor<{self.particle}>"
...@@ -62,11 +61,15 @@ class NeighborFor: ...@@ -62,11 +61,15 @@ class NeighborFor:
particle_shape = self.sim.particle_shape particle_shape = self.sim.particle_shape
nshapes = self.cell_lists.nshapes nshapes = self.cell_lists.nshapes
for shape in range(self.sim.max_shapes()): for shape in self.shapes:
for disp in For(self.sim, 0, self.cell_lists.nstencil): for disp in For(self.sim, -1, self.cell_lists.nstencil):
neigh_cell = particle_cell[self.particle] + stencil[disp] neigh_cell = \
Select(self.sim, disp < 0, 0, particle_cell[self.particle] + stencil[disp])
for _ in Filter(self.sim, ScalarOp.or_op(disp < 0,
ScalarOp.and_op(neigh_cell > 0,
neigh_cell < ncells))):
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) 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]): for cell_particle in For(self.sim, start, start + nshapes[neigh_cell][shape]):
particle_id = cell_particles[neigh_cell][cell_particle] particle_id = cell_particles[neigh_cell][cell_particle]
...@@ -81,26 +84,11 @@ class NeighborFor: ...@@ -81,26 +84,11 @@ class NeighborFor:
for _ in Filter(self.sim, condition): for _ in Filter(self.sim, condition):
yield Neighbor(self.sim, cell_particle, neigh_cell, particle_id, shape) 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]
if self.sim._compute_half:
shape_cond = particle_shape[inf_particle] > self.sim.get_shape_id(shape)
condition = ScalarOp.or_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: else:
neighborlists = self.neighbor_lists.neighborlists neighborlists = self.neighbor_lists.neighborlists
numneighs = self.neighbor_lists.numneighs numneighs = self.neighbor_lists.numneighs
for shape in range(self.sim.max_shapes()): for shape in self.shapes:
start = sum([numneighs[self.particle][s] for s in range(shape)], 0) start = sum([numneighs[self.particle][s] for s in range(shape)], 0)
for k in For(self.sim, start, start + numneighs[self.particle][shape]): for k in For(self.sim, start, start + numneighs[self.particle][shape]):
yield Neighbor(self.sim, k, None, neighborlists[self.particle][k], shape) yield Neighbor(self.sim, k, None, neighborlists[self.particle][k], shape)
...@@ -143,17 +131,18 @@ class InteractionData: ...@@ -143,17 +131,18 @@ class InteractionData:
class ParticleInteraction(Lowerable): class ParticleInteraction(Lowerable):
def __init__(self, sim, nbody, cutoff_radius, use_cell_lists=False): def __init__(self, sim, nbody, cutoff_radius, use_cell_lists=False, split_kernels=False):
super().__init__(sim) super().__init__(sim)
self.nbody = nbody self.nbody = nbody
self.cutoff_radius = cutoff_radius self.cutoff_radius = cutoff_radius
self.contact_threshold = 0.0 self.contact_threshold = 0.0
self.use_cell_lists = use_cell_lists self.use_cell_lists = use_cell_lists
self.ncases = NeighborFor.number_of_cases(sim, use_cell_lists) self.split_kernels = split_kernels
self.interactions_data = [] self.nkernels = sim.max_shapes() if split_kernels else 1
self.blocks = [Block(sim, []) for _ in range(sim.max_shapes() * self.ncases)] self.interactions_data = {}
self.blocks = [Block(sim, []) for _ in range(sim.max_shapes())]
self.apply_list = [set() for _ in range(self.nkernels)]
self.active_block = None self.active_block = None
self.apply_list = set()
def add_statement(self, stmt): def add_statement(self, stmt):
self.active_block.add_statement(stmt) self.active_block.add_statement(stmt)
...@@ -161,16 +150,16 @@ class ParticleInteraction(Lowerable): ...@@ -161,16 +150,16 @@ class ParticleInteraction(Lowerable):
def __iter__(self): def __iter__(self):
self.sim.add_statement(self) self.sim.add_statement(self)
self.sim.enter(self) self.sim.enter(self)
self.sim.use_apply_list(self.apply_list)
# Neighbors vary across iterations # Neighbors vary across iterations
for shape in range(self.sim.max_shapes()): for shape in range(self.sim.max_shapes()):
for case in range(self.ncases): apply_list_id = shape if self.split_kernels else 0
self.active_block = self.blocks[shape * self.ncases + case] self.sim.use_apply_list(self.apply_list[apply_list_id])
self.interactions_data.append(InteractionData(self.sim, shape)) self.active_block = self.blocks[shape]
yield self.interactions_data[-1] self.interactions_data[shape] = InteractionData(self.sim, shape)
yield self.interactions_data[shape]
self.sim.release_apply_list()
self.sim.release_apply_list()
self.sim.leave() self.sim.leave()
self.active_block = None self.active_block = None
...@@ -181,88 +170,91 @@ class ParticleInteraction(Lowerable): ...@@ -181,88 +170,91 @@ class ParticleInteraction(Lowerable):
cell_lists = self.sim.cell_lists cell_lists = self.sim.cell_lists
neighbor_lists = None if self.use_cell_lists else self.sim.neighbor_lists neighbor_lists = None if self.use_cell_lists else self.sim.neighbor_lists
for i in ParticleFor(self.sim): for kernel in range(self.nkernels):
for _ in Filter(self.sim, ScalarOp.cmp(self.sim.particle_flags[i] & Flags.Fixed, 0)): for i in ParticleFor(self.sim):
for app in self.apply_list: for _ in Filter(self.sim, ScalarOp.cmp(self.sim.particle_flags[i] & Flags.Fixed, 0)):
app.add_reduction_variable() for app in self.apply_list[kernel]:
app.add_reduction_variable()
interaction = 0
for neigh in NeighborFor(self.sim, i, cell_lists, neighbor_lists): shapes = [kernel] if self.split_kernels else None
interaction_data = self.interactions_data[interaction] interaction = kernel
shape = interaction_data.shape() for neigh in NeighborFor(self.sim, i, cell_lists, neighbor_lists, shapes):
shape_id = self.sim.get_shape_id(shape) interaction_data = self.interactions_data[interaction]
j = neigh.particle_index() shape = interaction_data.shape()
shape_id = self.sim.get_shape_id(shape)
if shape_id == Shapes.PointMass: j = neigh.particle_index()
delta = position[i] - position[j]
squared_distance = delta.x() * delta.x() + \ if shape_id == Shapes.PointMass:
delta.y() * delta.y() + \ delta = position[i] - position[j]
delta.z() * delta.z() squared_distance = delta.x() * delta.x() + \
separation_dist = self.cutoff_radius * self.cutoff_radius delta.y() * delta.y() + \
cutoff_condition = squared_distance < separation_dist delta.z() * delta.z()
distance = Sqrt(self.sim, squared_distance) separation_dist = self.cutoff_radius * self.cutoff_radius
penetration_depth = None cutoff_condition = squared_distance < separation_dist
contact_normal = None distance = Sqrt(self.sim, squared_distance)
contact_point = None penetration_depth = None
contact_normal = None
elif shape_id == Shapes.Sphere: contact_point = None
radius = self.sim.property('radius')
delta = position[i] - position[j] elif shape_id == Shapes.Sphere:
squared_distance = delta.x() * delta.x() + \ radius = self.sim.property('radius')
delta.y() * delta.y() + \ delta = position[i] - position[j]
delta.z() * delta.z() squared_distance = delta.x() * delta.x() + \
separation_dist = radius[i] + radius[j] + self.contact_threshold delta.y() * delta.y() + \
cutoff_condition = squared_distance < separation_dist * separation_dist delta.z() * delta.z()
distance = Sqrt(self.sim, squared_distance) separation_dist = radius[i] + radius[j] + self.contact_threshold
penetration_depth = distance - radius[i] - radius[j] cutoff_condition = squared_distance < separation_dist * separation_dist
contact_normal = delta * (1.0 / distance) distance = Sqrt(self.sim, squared_distance)
k = radius[j] + 0.5 * penetration_depth penetration_depth = distance - radius[i] - radius[j]
contact_point = position[j] + contact_normal * k contact_normal = delta * (1.0 / distance)
k = radius[j] + 0.5 * penetration_depth
elif shape_id == Shapes.Halfspace: contact_point = position[j] + contact_normal * k
radius = self.sim.property('radius')
normal = self.sim.property('normal') elif shape_id == Shapes.Halfspace:
radius = self.sim.property('radius')
d = normal[j][0] * position[j][0] + \ normal = self.sim.property('normal')
normal[j][1] * position[j][1] + \
normal[j][2] * position[j][2] d = normal[j][0] * position[j][0] + \
normal[j][1] * position[j][1] + \
k = normal[j][0] * position[i][0] + \ normal[j][2] * position[j][2]
normal[j][1] * position[i][1] + \
normal[j][2] * position[i][2] k = normal[j][0] * position[i][0] + \
normal[j][1] * position[i][1] + \
penetration_depth = k - radius[i] - d normal[j][2] * position[i][2]
cutoff_condition = penetration_depth < self.contact_threshold
tmp = radius[i] + penetration_depth penetration_depth = k - radius[i] - d
contact_normal = normal[j] cutoff_condition = penetration_depth < self.contact_threshold
contact_point = position[i] - Vector(self.sim, [tmp, tmp, tmp]) * normal[j] tmp = radius[i] + penetration_depth
contact_normal = normal[j]
else: contact_point = position[i] - Vector(self.sim, [tmp, tmp, tmp]) * normal[j]
raise Exception("Invalid shape id.")
else:
interaction_data.i().assign(i) raise Exception("Invalid shape id.")
interaction_data.j().assign(j)
interaction_data.delta().assign(delta) interaction_data.i().assign(i)
interaction_data.squared_distance().assign(squared_distance) interaction_data.j().assign(j)
interaction_data.penetration_depth().assign(penetration_depth) interaction_data.delta().assign(delta)
interaction_data.contact_point().assign(contact_point) interaction_data.squared_distance().assign(squared_distance)
interaction_data.contact_normal().assign(contact_normal) interaction_data.penetration_depth().assign(penetration_depth)
self.sim.add_statement(Filter(self.sim, cutoff_condition, self.blocks[interaction])) interaction_data.contact_point().assign(contact_point)
interaction += 1 interaction_data.contact_normal().assign(contact_normal)
self.sim.add_statement(
prop_reductions = {} Filter(self.sim, cutoff_condition, self.blocks[shape]))
for app in self.apply_list: interaction += 1
prop = app.prop()
reduction = app.reduction_variable() prop_reductions = {}
for app in self.apply_list[kernel]:
if prop not in prop_reductions: prop = app.prop()
prop_reductions[prop] = reduction reduction = app.reduction_variable()
else: if prop not in prop_reductions:
prop_reductions[prop] = prop_reductions[prop] + reduction prop_reductions[prop] = reduction
for prop, reduction in prop_reductions.items(): else:
Assign(self.sim, prop[i], prop[i] + reduction) prop_reductions[prop] = prop_reductions[prop] + reduction
for prop, reduction in prop_reductions.items():
Assign(self.sim, prop[i], prop[i] + reduction)
else: else:
raise Exception("Interactions among more than two particles are currently not supported.") raise Exception("Interactions among more than two particles are currently not supported.")
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment