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

Fix half-neighbor lists

parent 7b78e4a8
No related branches found
No related tags found
1 merge request!1Implement DEM and many other features
......@@ -5,7 +5,7 @@ import sys
def lj(i, j):
sr2 = 1.0 / squared_distance(i, j)
sr6 = sr2 * sr2 * sr2 * sigma6[i, j]
apply(force, delta(i, j) * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j])
apply(force, delta(i, j) * (48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[i, j]))
def euler(i):
linear_velocity[i] += dt * force[i] / mass[i]
......@@ -36,6 +36,7 @@ psim.add_feature_property('type', 'sigma6', pairs.real(), [epsilon for i in rang
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.compute_half()
psim.build_neighbor_lists(cutoff_radius + skin)
psim.vtk_output(f"output/lj_{target}")
psim.compute(lj, cutoff_radius)
......
......@@ -44,11 +44,11 @@ class Apply(Lowerable):
@pairs_inline
def lower(self):
Assign(self.sim, self._reduction_variable, self._reduction_variable + self._expr)
if self.sim._compute_half:
for _ in Filter(self.sim,
ScalarOp.and_op(self._j < self.sim.nlocal,
ScalarOp.cmp(self.sim.particle_flags[self._j] & Flags.Fixed, 0))):
Assign(self.sim, self._prop[self._j], self._prop[self._j] - self._expr)
Assign(self.sim, self._reduction_variable, self._reduction_variable + self._expr)
......@@ -72,10 +72,8 @@ class NeighborFor:
particle_id = cell_particles[neigh_cell][cell_particle]
if self.sim._compute_half:
shape_cond = particle_shape[particle_id] > shape
condition = ScalarOp.or_op(shape_cond,
ScalarOp.and_op(ScalarOp.not_op(shape_cond),
self.particle < particle_id))
shape_cond = particle_shape[particle_id] > self.sim.get_shape_id(shape)
condition = ScalarOp.or_op(shape_cond, self.particle < particle_id)
else:
condition = ScalarOp.neq(particle_id, self.particle)
......@@ -89,10 +87,8 @@ class NeighborFor:
inf_particle = cell_particles[0][inf_id]
if self.sim._compute_half:
shape_cond = particle_shape[inf_particle] > shape
condition = ScalarOp.or_op(shape_cond,
ScalarOp.and_op(ScalarOp.not_op(shape_cond),
self.particle < inf_particle))
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)
......@@ -266,7 +262,11 @@ class ParticleInteraction(Lowerable):
prop_reductions[prop] = prop_reductions[prop] + reduction
for prop, reduction in prop_reductions.items():
Assign(self.sim, prop[i], reduction)
if self.sim._compute_half:
Assign(self.sim, prop[i], prop[i] + reduction)
else:
Assign(self.sim, prop[i], reduction)
else:
raise Exception("Interactions among more than two particles are currently not supported.")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment