Merge pull request #189 from prabhuramachandran/spatial-ordering-option

Fix bug in spatial ordering and add option for it.

pysph/base/nnps_base.pyx

@@ -1397,5 +1397,6 @@ cdef class NNPS(NNPSBase):
         self.get_spatially_ordered_indices(pa_index, indices)
         cdef BaseArray arr
 
-        for arr in pa.properties.values():
-            arr.c_align_array(indices)
+        for name, arr in pa.properties.items():
+            stride = pa.stride.get(name, 1)
+            arr.c_align_array(indices, stride)

pysph/solver/application.py

@@ -376,6 +376,12 @@ def _setup_argparse(self):
             type=int,
             help="Printing frequency for the output")
 
+        parser.add_argument(
+            '--reorder-freq', action="store", dest="reorder_freq",
+            default=None, type=int,
+            help="Frequency between spatially reordering particles."
+        )
+
         # --detailed-output.
         parser.add_argument(
             "--detailed-output",
@@ -1123,6 +1129,12 @@ def _configure_solver(self):
         # disable_output
         solver.set_disable_output(options.disable_output)
 
+        if options.reorder_freq is None:
+            if options.with_opencl:
+                solver.set_reorder_freq(50)
+        else:
+            solver.set_reorder_freq(options.reorder_freq)
+
         # output print frequency
         if options.freq is not None:
             solver.set_print_freq(options.freq)

pysph/solver/solver.py

@@ -69,6 +69,10 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         fixed_h : bint
             Flag for constant smoothing lengths `h`
 
+        reorder_freq : int
+            The number of iterations after which particles should
+            be re-ordered.  If zero, do not do this.
+
         Example
         -------
 
@@ -91,6 +95,7 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         self.particles = None
 
         self.acceleration_evals = None
+        self.nnps = None
 
         # solver time and iteration count
         self.t = 0
@@ -164,6 +169,8 @@ def __init__(self, dim=2, integrator=None, kernel=None,
         # flag for constant smoothing lengths
         self.fixed_h = fixed_h
 
+        self.reorder_freq = 0
+
         # Set all extra keyword arguments
         for attr, value in kwargs.items():
             if hasattr(self, attr):
@@ -203,6 +210,7 @@ def setup(self, particles, equations, nnps, kernel=None, fixed_h=False):
         sph_compiler.compile()
 
         # Set the nnps for all concerned objects.
+        self.nnps = nnps
         for ae in self.acceleration_evals:
             ae.set_nnps(nnps)
         self.integrator.set_nnps(nnps)
@@ -283,6 +291,14 @@ def append_particle_arrrays(self, arrays):
 
         self.setup(self.particles)
 
+    def reorder_particles(self):
+        """Re-order particles so as to coalesce memory access.
+        """
+        for i in range(len(self.particles)):
+            self.nnps.spatially_order_particles(i)
+        # We must update after the reorder.
+        self.nnps.update()
+
     def set_adaptive_timestep(self, value):
         """Set it to True to use adaptive timestepping based on
         cfl, viscous and force factor.
@@ -395,6 +411,11 @@ def set_command_handler(self, callable, command_interval=1):
     def set_parallel_manager(self, pm):
         self.pm = pm
 
+    def set_reorder_freq(self, freq):
+        """Set the reorder frequency in number of iterations.
+        """
+        self.reorder_freq = freq
+
     def barrier(self):
         if self.comm:
             self.comm.barrier()
@@ -422,6 +443,10 @@ def solve(self, show_progress=True):
         self.dump_output()
         self.barrier()  # everybody waits for this to complete
 
+        reorder_freq = self.reorder_freq
+        if reorder_freq > 0:
+            self.reorder_particles()
+
         # Compute the accelerations once for the predictor corrector
         # integrator to work correctly at the first time step.
         self.integrator.initial_acceleration(self.t, self.dt)
@@ -468,6 +493,9 @@ def solve(self, show_progress=True):
             # update the time for all arrays
             self.update_particle_time()
 
+            if reorder_freq > 0 and (self.count % reorder_freq == 0):
+                self.reorder_particles()
+
             if self.execute_commands is not None:
                 if self.count % self.command_interval == 0:
                     self.execute_commands(self)