From ffa2436b478796506e3ac6eae39a95d283b30643 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Thu, 9 Nov 2023 03:12:42 +0100
Subject: [PATCH] Update LJ example
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
examples/dem.py | 8 +++++++-
examples/lj.py | 13 ++++++-------
src/pairs/__init__.py | 7 +++++--
src/pairs/mapping/keywords.py | 6 ++++++
src/pairs/sim/cell_lists.py | 8 +++++++-
src/pairs/sim/comm.py | 4 ++--
src/pairs/sim/interaction.py | 25 ++++++++++++++++++++-----
src/pairs/sim/shapes.py | 2 +-
src/pairs/sim/simulation.py | 8 ++++++--
9 files changed, 60 insertions(+), 21 deletions(-)
diff --git a/examples/dem.py b/examples/dem.py
index 557440c..0ab64a3 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -131,7 +131,13 @@ lnDryResCoeff = math.log(restitutionCoefficient);
frictionStatic = 0.0
frictionDynamic = frictionCoefficient
-psim = pairs.simulation("dem", timesteps=timeSteps, double_prec=True, use_contact_history=True)
+psim = pairs.simulation(
+ "dem",
+ [pairs.sphere(), pairs.halfspace()],
+ timesteps=timeSteps,
+ double_prec=True,
+ use_contact_history=True)
+
#psim = pairs.simulation("dem", debug=True, timesteps=timeSteps)
psim.add_position('position')
psim.add_property('mass', pairs.real(), 1.0)
diff --git a/examples/lj.py b/examples/lj.py
index fafdfd4..be5883c 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -26,17 +26,16 @@ sigma = 1.0
epsilon = 1.0
sigma6 = sigma ** 6
-psim = pairs.simulation("lj", debug=True, timesteps=200)
+psim = pairs.simulation("lj", [pairs.point_mass()], timesteps=200, double_prec=True)
psim.add_position('position')
-psim.add_property('mass', pairs.double(), 1.0)
+psim.add_property('mass', pairs.real(), 1.0)
psim.add_property('linear_velocity', pairs.vector())
psim.add_property('force', pairs.vector(), volatile=True)
psim.add_feature('type', ntypes)
-psim.add_feature_property('type', 'epsilon', pairs.double(), [sigma for i in range(ntypes * ntypes)])
-psim.add_feature_property('type', 'sigma6', pairs.double(), [epsilon for i in range(ntypes * ntypes)])
-#psim.set_domain(0.0, 0.0, 0.0, 6.7184, 6.7184, 6.7184) # for 4x4x4 setup
-psim.set_domain(0.0, 0.0, 0.0, 53.747078, 53.747078, 53.747078)
-psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'linear_velocity', 'flags'], pairs.sphere())
+psim.add_feature_property('type', 'epsilon', pairs.real(), [sigma for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'sigma6', pairs.real(), [epsilon for i in range(ntypes * ntypes)])
+psim.set_domain([0.0, 0.0, 0.0, 53.747078, 53.747078, 53.747078])
+psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'linear_velocity', 'flags'], pairs.point_mass())
psim.reneighbor_every(20)
psim.build_neighbor_lists(cutoff_radius + skin)
psim.vtk_output(f"output/lj_{target}")
diff --git a/src/pairs/__init__.py b/src/pairs/__init__.py
index b73e6e9..7ab7b1a 100644
--- a/src/pairs/__init__.py
+++ b/src/pairs/__init__.py
@@ -5,8 +5,8 @@ from pairs.sim.shapes import Shapes
from pairs.sim.simulation import Simulation
-def simulation(ref, dims=3, timesteps=100, double_prec=False, use_contact_history=False, debug=False):
- return Simulation(CGen(ref, debug), dims, timesteps, double_prec, use_contact_history)
+def simulation(ref, shapes, dims=3, timesteps=100, double_prec=False, use_contact_history=False, debug=False):
+ return Simulation(CGen(ref, debug), shapes, dims, timesteps, double_prec, use_contact_history)
def target_cpu():
return Target(Target.Backend_CPP, Target.Feature_CPU)
@@ -35,6 +35,9 @@ def matrix():
def quaternion():
return Types.Quaternion
+def point_mass():
+ return Shapes.PointMass
+
def sphere():
return Shapes.Sphere
diff --git a/src/pairs/mapping/keywords.py b/src/pairs/mapping/keywords.py
index 3754454..2d7917f 100644
--- a/src/pairs/mapping/keywords.py
+++ b/src/pairs/mapping/keywords.py
@@ -29,6 +29,12 @@ class Keywords:
method = self.get_method(f"keyword_{keyword}")
return method is not None
+ def keyword_is_point_mass(self, args):
+ assert len(args) == 1, "is_point_mass() keyword requires one parameter."
+ particle_id = args[0]
+ assert particle_id.type() == Types.Int32, "Particle ID must be an integer."
+ return ScalarOp.cmp(self.sim.particle_shape[particle_id], Shapes.PointMass)
+
def keyword_is_sphere(self, args):
assert len(args) == 1, "is_sphere() keyword requires one parameter."
particle_id = args[0]
diff --git a/src/pairs/sim/cell_lists.py b/src/pairs/sim/cell_lists.py
index 955f622..697ab87 100644
--- a/src/pairs/sim/cell_lists.py
+++ b/src/pairs/sim/cell_lists.py
@@ -29,6 +29,7 @@ class CellLists:
self.ncells_capacity = self.sim.add_var('ncells_capacity', Types.Int32, 100)
self.cell_capacity = self.sim.add_var('cell_capacity', Types.Int32, 20)
self.dim_ncells = self.sim.add_static_array('dim_cells', self.sim.ndims(), Types.Int32)
+ self.shapes_buffer = self.sim.add_static_array('shapes_buffer', self.sim.max_shapes(), Types.Int32)
self.cell_particles = self.sim.add_array('cell_particles', [self.ncells_capacity, self.cell_capacity], Types.Int32)
self.cell_sizes = self.sim.add_array('cell_sizes', self.ncells_capacity, Types.Int32)
self.nshapes = self.sim.add_array('nshapes', [self.ncells_capacity, self.sim.max_shapes()], Types.Int32)
@@ -116,19 +117,24 @@ class PartitionCellLists(Lowerable):
def lower(self):
self.sim.module_name("partition_cell_lists")
cell_particles = self.cell_lists.cell_particles
+ shapes_buffer = self.cell_lists.shapes_buffer
+
+ for s in range(self.sim.max_shapes()):
+ Assign(self.sim, shapes_buffer[s], self.sim.get_shape_id(s))
for cell in For(self.sim, 0, self.cell_lists.ncells):
start = self.sim.add_temp_var(0)
end = self.sim.add_temp_var(0)
for shape in For(self.sim, 0, self.sim.max_shapes()):
+ shape_id = shapes_buffer[shape]
shape_start = self.sim.add_temp_var(start)
Assign(self.sim, end, self.cell_lists.cell_sizes[cell] - 1)
for _ in While(self.sim, start <= end):
particle = cell_particles[cell][start]
- for unmatch in Branch(self.sim, ScalarOp.neq(self.sim.particle_shape[particle], shape)):
+ for unmatch in Branch(self.sim, ScalarOp.neq(self.sim.particle_shape[particle], shape_id)):
if unmatch:
for _ in Filter(self.sim, ScalarOp.neq(start, end)):
Assign(self.sim, cell_particles[cell][start], cell_particles[cell][end])
diff --git a/src/pairs/sim/comm.py b/src/pairs/sim/comm.py
index af60379..7e114ef 100644
--- a/src/pairs/sim/comm.py
+++ b/src/pairs/sim/comm.py
@@ -38,7 +38,7 @@ class Comm:
def synchronize(self):
# Every property that is not 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 ['position', 'linear_velocity', 'angular_velocity']]
+ prop_list = [self.sim.property(p) for p in ['position', 'linear_velocity', 'angular_velocity'] if self.sim.property(p) is not None]
for step in range(self.dom_part.number_of_steps()):
PackGhostParticles(self, step, prop_list)
CommunicateData(self, step, prop_list)
@@ -49,7 +49,7 @@ class Comm:
# 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 ['uid', 'mass', 'radius', 'position', 'linear_velocity', 'angular_velocity', 'shape']]
+ prop_list = [self.sim.property(p) for p in ['uid', 'mass', 'radius', 'position', 'linear_velocity', 'angular_velocity', 'shape'] if self.sim.property(p) is not None]
Assign(self.sim, self.nsend_all, 0)
Assign(self.sim, self.sim.nghost, 0)
diff --git a/src/pairs/sim/interaction.py b/src/pairs/sim/interaction.py
index f76fdbf..02bf67f 100644
--- a/src/pairs/sim/interaction.py
+++ b/src/pairs/sim/interaction.py
@@ -158,8 +158,6 @@ class ParticleInteraction(Lowerable):
def lower(self):
if self.nbody == 2:
position = self.sim.position()
- radius = self.sim.property('radius')
- normal = self.sim.property('normal')
cell_lists = self.sim.cell_lists
neighbor_lists = None if self.use_cell_lists else self.sim.neighbor_lists
@@ -169,9 +167,23 @@ class ParticleInteraction(Lowerable):
for neigh in NeighborFor(self.sim, i, cell_lists, neighbor_lists):
interaction_data = self.interactions_data[interaction]
shape = interaction_data.shape()
+ shape_id = self.sim.get_shape_id(shape)
j = neigh.particle_index()
- if shape == Shapes.Sphere:
+ if shape_id == Shapes.PointMass:
+ delta = position[i] - position[j]
+ squared_distance = delta.x() * delta.x() + \
+ delta.y() * delta.y() + \
+ delta.z() * delta.z()
+ separation_dist = self.cutoff_radius
+ cutoff_condition = squared_distance < self.cutoff_radius
+ distance = Sqrt(self.sim, squared_distance)
+ penetration_depth = None
+ contact_normal = None
+ contact_point = None
+
+ elif shape_id == Shapes.Sphere:
+ radius = self.sim.property('radius')
delta = position[i] - position[j]
squared_distance = delta.x() * delta.x() + \
delta.y() * delta.y() + \
@@ -184,7 +196,10 @@ class ParticleInteraction(Lowerable):
k = radius[j] + 0.5 * penetration_depth
contact_point = position[j] + contact_normal * k
- elif shape == Shapes.Halfspace:
+ elif shape_id == Shapes.Halfspace:
+ radius = self.sim.property('radius')
+ normal = self.sim.property('normal')
+
d = normal[j][0] * position[j][0] + \
normal[j][1] * position[j][1] + \
normal[j][2] * position[j][2]
@@ -200,7 +215,7 @@ class ParticleInteraction(Lowerable):
contact_point = position[i] - Vector(self.sim, [tmp, tmp, tmp]) * normal[j]
else:
- raise Exception("Invalid shape!")
+ raise Exception("Invalid shape id.")
interaction_data.i().assign(i)
interaction_data.j().assign(j)
diff --git a/src/pairs/sim/shapes.py b/src/pairs/sim/shapes.py
index dbd37db..3f768c6 100644
--- a/src/pairs/sim/shapes.py
+++ b/src/pairs/sim/shapes.py
@@ -1,4 +1,4 @@
class Shapes:
Sphere = 0
Halfspace = 1
- PointMasses = 2
+ PointMass = 2
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index d55722d..302868c 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -32,7 +32,7 @@ from pairs.transformations import Transformations
class Simulation:
- def __init__(self, code_gen, dims=3, timesteps=100, double_prec=False, use_contact_history=False):
+ def __init__(self, code_gen, shapes, dims=3, timesteps=100, double_prec=False, use_contact_history=False):
self.code_gen = code_gen
self.code_gen.assign_simulation(self)
self.position_prop = None
@@ -79,12 +79,16 @@ class Simulation:
self._pbc = [True for _ in range(dims)]
self._use_contact_history = use_contact_history
self._contact_history = ContactHistory(self) if use_contact_history else None
+ self._shapes = shapes
def use_double_precision(self):
return self._double_prec
+ def get_shape_id(self, shape):
+ return self._shapes[shape]
+
def max_shapes(self):
- return 2
+ return len(self._shapes)
def add_module(self, module):
assert isinstance(module, Module), "add_module(): Given parameter is not of type Module!"
--
GitLab