diff --git a/examples/dem.py b/examples/dem.py
index 33a4e366756b7419884a95cd834b856b8c9b3900..5e7a6a02852183d450243db2381c6d6105b2d583 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -146,6 +146,7 @@ psim.add_contact_property('tangential_spring_displacement', pairs.vector(), [0.0
psim.add_contact_property('impact_velocity_magnitude', pairs.double(), 0.0)
psim.set_domain([0.0, 0.0, 0.0, domainSize_SI[0], domainSize_SI[1], domainSize_SI[2]])
+psim.pbc([True, True, False])
psim.read_particle_data(
"data/spheres.input",
['type', 'mass', 'radius', 'position', 'linear_velocity', 'flags'],
@@ -163,12 +164,12 @@ psim.read_particle_data(
psim.build_neighbor_lists(linkedCellWidth + skin)
psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
+#psim.compute(gravity, symbols={'densityParticle_SI': densityParticle_SI,
+# 'densityFluid_SI': densityFluid_SI,
+# 'gravity_SI': gravity_SI,
+# 'pi': math.pi })
psim.compute(linear_spring_dashpot, linkedCellWidth + skin, symbols={'dt': dt_SI})
psim.compute(euler, symbols={'dt': dt_SI})
-psim.compute(gravity, symbols={'densityParticle_SI': densityParticle_SI,
- 'densityFluid_SI': densityFluid_SI,
- 'gravity_SI': gravity_SI,
- 'pi': math.pi })
if target == 'gpu':
psim.target(pairs.target_gpu())
diff --git a/examples/lj.py b/examples/lj.py
index c5989fe8b7b37b73912b7f8b0d331aef02335e48..fafdfd4e7edd8331f0b2e53a42606ae8edf6f936 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -37,6 +37,7 @@ psim.add_feature_property('type', 'sigma6', pairs.double(), [epsilon for i in ra
#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.reneighbor_every(20)
psim.build_neighbor_lists(cutoff_radius + skin)
psim.vtk_output(f"output/lj_{target}")
psim.compute(lj, cutoff_radius)
diff --git a/src/pairs/mapping/funcs.py b/src/pairs/mapping/funcs.py
index 2e753d70dfbca91805a4c4e9afd4eec292f7e443..c57a209bbe70718a9bb7115060662363476748da 100644
--- a/src/pairs/mapping/funcs.py
+++ b/src/pairs/mapping/funcs.py
@@ -101,11 +101,13 @@ class BuildParticleIR(ast.NodeVisitor):
def visit_AugAssign(self, node):
lhs = self.visit(node.target)
rhs = self.visit(node.value)
+ op_class = VectorOp if Types.Vector in [lhs.type(), rhs.type()] else ScalarOp
+ bin_op = op_class(self.sim, lhs, rhs, BuildParticleIR.get_binary_op(node.op))
if isinstance(lhs, UndefinedSymbol):
- self.add_symbols({lhs.symbol_id: rhs})
+ self.add_symbols({lhs.symbol_id: bin_op})
else:
- Assign(self.sim, lhs, lhs + rhs)
+ Assign(self.sim, lhs, bin_op)
def visit_BinOp(self, node):
#print(ast.dump(node))
diff --git a/src/pairs/sim/domain_partitioning.py b/src/pairs/sim/domain_partitioning.py
index 86aa15b7f5dbb77726f3c3c4367746c852b5555b..08ce9a48d41f25cf8524f1a5665c3ef638b416ab 100644
--- a/src/pairs/sim/domain_partitioning.py
+++ b/src/pairs/sim/domain_partitioning.py
@@ -23,6 +23,9 @@ class DimensionRanges:
return [step * 2 + 0, step * 2 + 1]
def ghost_particles(self, step, position, offset=0.0):
+ if self.sim._pbc[step] is False:
+ return
+
# Particles with one of the following flags are ignored
flags_to_exclude = (Flags.Infinite | Flags.Fixed | Flags.Global)
diff --git a/src/pairs/sim/pbc.py b/src/pairs/sim/pbc.py
index 277c18f3f373e65861e282c7979d87fe894c25ac..af5fec1516b8aa485adbbcd8f578cd93bbb6ccea 100644
--- a/src/pairs/sim/pbc.py
+++ b/src/pairs/sim/pbc.py
@@ -23,8 +23,9 @@ class EnforcePBC(Lowerable):
for i in ParticleFor(sim):
# TODO: VecFilter?
for d in range(0, ndims):
- for _ in Filter(sim, positions[i][d] < grid.min(d)):
- Assign(sim, positions[i][d], positions[i][d] + grid.length(d))
+ if sim._pbc[d] is True:
+ for _ in Filter(sim, positions[i][d] < grid.min(d)):
+ Assign(sim, positions[i][d], positions[i][d] + grid.length(d))
- for _ in Filter(sim, positions[i][d] > grid.max(d)):
- Assign(sim, positions[i][d], positions[i][d] - grid.length(d))
+ for _ in Filter(sim, positions[i][d] > grid.max(d)):
+ Assign(sim, positions[i][d], positions[i][d] - grid.length(d))
diff --git a/src/pairs/sim/shapes.py b/src/pairs/sim/shapes.py
index 8359da7f2b6a1f7fd138d100a3eb77035d70e4ed..dbd37db4107e59a53577413901c1a33a33467b14 100644
--- a/src/pairs/sim/shapes.py
+++ b/src/pairs/sim/shapes.py
@@ -1,3 +1,4 @@
class Shapes:
Sphere = 0
Halfspace = 1
+ PointMasses = 2
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 6315036dd311555cbb5ddcc00f44ebc3059c7fa5..8fd8b85e87cef5a8882773f46a29ec5f7e697eb3 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -68,10 +68,12 @@ class Simulation:
self.ntimesteps = timesteps
self.expr_id = 0
self.iter_id = 0
+ self.reneighbor_frequency = 1
self.vtk_file = None
self.vtk_frequency = 0
self._target = None
self._dom_part = DimensionRanges(self)
+ self._pbc = [True for _ in range(dims)]
def max_shapes(self):
return 2
@@ -107,6 +109,10 @@ class Simulation:
def ndims(self):
return self.dims
+ def pbc(self, pbc_config):
+ assert len(pbc_config) == self.dims, "PBC must be specified for each dimension."
+ self._pbc = pbc_config
+
def add_property(self, prop_name, prop_type, value=0.0, volatile=False):
assert self.property(prop_name) is None, f"Property already defined: {prop_name}"
return self.properties.add(prop_name, prop_type, value, volatile)
@@ -173,6 +179,9 @@ class Simulation:
self.grid = Grid3D(self, grid[0], grid[1], grid[2], grid[3], grid[4], grid[5])
self.setups.add_statement(InitializeDomain(self))
+ def reneighbor_every(self, frequency):
+ self.reneighbor_frequency = frequency
+
def create_particle_lattice(self, grid, spacing, props={}):
self.setups.add_statement(ParticleLattice(self, grid, spacing, props, self.position()))
@@ -267,16 +276,16 @@ class Simulation:
comm = Comm(self, self._dom_part)
timestep_procedures = [
- (comm.exchange(), 20),
- (comm.borders(), comm.synchronize(), 20),
- (BuildCellLists(self, self.cell_lists), 20),
- (PartitionCellLists(self, self.cell_lists), 20),
- (BuildNeighborLists(self, self.neighbor_lists), 20),
+ (comm.exchange(), self.reneighbor_frequency),
+ (comm.borders(), comm.synchronize(), self.reneighbor_frequency),
+ (BuildCellLists(self, self.cell_lists), self.reneighbor_frequency),
+ (PartitionCellLists(self, self.cell_lists), self.reneighbor_frequency),
+ (BuildNeighborLists(self, self.neighbor_lists), self.reneighbor_frequency),
]
if not self.contact_properties.empty():
contact_history = ContactHistory(self.cell_lists)
- timestep_procedures.append((BuildContactHistory(self, contact_history), 20))
+ timestep_procedures.append((BuildContactHistory(self, contact_history), self.reneighbor_frequency))
timestep_procedures += [
ResetVolatileProperties(self),