Add options to change PBC and reneighbor frequency

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>

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())

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)

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))

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)
 

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)) 

src/pairs/sim/shapes.py

 class Shapes:
     Sphere      =   0
     Halfspace   =   1
+    PointMasses =   2

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),