Merge pull request #180 from prabhuramachandran/auto-update-domain

API: Integrators explicitly call update_domain.

docs/source/design/equations.rst

@@ -644,14 +644,22 @@ test suite to illustrate this:
 
             self.compute_accelerations(1, update_nnps=False)
             self.stage2()
+            self.update_domain()
+            self.do_post_stage(dt, 2)
 
 Note that the ``compute_accelerations`` method takes two arguments, the
 ``index`` (which defaults to zero) and ``update_nnps`` which defaults to
 ``True``. A simple integrator with a single RHS would simply call
 ``self.compute_accelerations()``. However, in the above, the first set of
 equations is called first, and then for the second stage the second set of
 equations is evaluated but without updating the NNPS (handy if the particles
-do not move in stage1).
+do not move in stage1). Note the call ``self.update_domain()`` after the
+second stage, this sets up any ghost particles for periodicity when particles
+have been moved, it also updates the neighbor finder to use an appropriate
+neighbor length based on the current smoothing length. If you do not need to
+do this for your particular integrator you may choose not to add this. In the
+above case, the domain is not updated after the first stage as the particles
+have not moved.
 
 The above illustrates how one can create more complex integrators that employ
 different accelerations in each stage.

pysph/base/gpu_domain_manager.py

@@ -32,7 +32,7 @@ def __init__(self, xmin=-1000., xmax=1000., ymin=0.,
 
         self.ghosts = None
 
-    def update(self, *args, **kwargs):
+    def update(self):
         """General method that is called before NNPS can bin particles.
 
         This method is responsible for the computation of cell sizes

pysph/base/gpu_nnps_base.pyx

@@ -249,7 +249,7 @@ cdef class GPUNNPS(NNPSBase):
             for cache in self.cache:
                 cache.update()
 
-    def update_domain(self, *args, **kwargs):
+    def update_domain(self):
         self.domain.update()
 
     cdef _compute_bounds(self):

pysph/base/nnps_base.pyx

@@ -439,7 +439,7 @@ cdef class CPUDomainManager(DomainManagerBase):
         self.ghosts = None
 
     #### Public protocol ################################################
-    def update(self, *args, **kwargs):
+    def update(self):
         """General method that is called before NNPS can bin particles.
 
         This method is responsible for the computation of cell sizes
@@ -488,7 +488,7 @@ cdef class CPUDomainManager(DomainManagerBase):
         """
         # minimum and maximum values of the domain
         cdef double xmin = self.xmin, xmax = self.xmax
-        cdef double ymin = self.ymin, ymax = self.ymax,
+        cdef double ymin = self.ymin, ymax = self.ymax
         cdef double zmin = self.zmin, zmax = self.zmax
 
         # translations along each coordinate direction
@@ -1238,7 +1238,7 @@ cdef class NNPS(NNPSBase):
     def set_in_parallel(self, bint in_parallel):
         self.domain.manager.in_parallel = in_parallel
 
-    def update_domain(self, *args, **kwargs):
+    def update_domain(self):
         self.domain.update()
 
     cpdef update(self):

pysph/solver/application.py

@@ -1080,11 +1080,6 @@ def _configure_solver(self):
             self.nnps = nnps
 
         nnps = self.nnps
-        # once the NNPS has been set-up, we set the default Solver
-        # post-stage callback to the DomainManager.setup_domain
-        # method. This method is responsible to computing the new cell
-        # size and doing any periodicity checks if needed.
-        solver.add_post_stage_callback(nnps.update_domain)
 
         # inform NNPS if it's working in parallel
         if self.num_procs > 1:
@@ -1463,7 +1458,6 @@ def run(self, argv=None):
         if self.options.with_opencl and self.options.profile:
             from compyle.opencl import print_profile
             print_profile()
-            print_mem_usage()
         self._write_info(
             self.info_filename, completed=True, cpu_time=run_duration)
 

pysph/sph/integrator.py

@@ -160,8 +160,9 @@ def one_timestep(self, t, dt):
                 - self.stage1(), self.stage2() etc. depending on the number of
                   stages available.
 
-                - self.compute_accelerations(t, dt)
+                - self.compute_accelerations(index=0, update_nnps=True)
                 - self.do_post_stage(stage_dt, stage_count_from_1)
+                - self.update_domain()
 
         Please see any of the concrete implementations of the Integrator class
         to study.  By default the Integrator implements a
@@ -172,6 +173,7 @@ def one_timestep(self, t, dt):
 
         # Predict
         self.stage1()
+        self.update_domain()
 
         # Call any post-stage functions.
         self.do_post_stage(0.5*dt, 1)
@@ -180,6 +182,7 @@ def one_timestep(self, t, dt):
 
         # Correct
         self.stage2()
+        self.update_domain()
 
         # Call any post-stage functions.
         self.do_post_stage(dt, 2)
@@ -236,6 +239,21 @@ def initial_acceleration(self, t, dt):
         """
         self.acceleration_evals[0].compute(t, dt)
 
+    def update_domain(self):
+        """Update the domain of the simulation.
+
+        This is to be called when particles move so the ghost particles
+        (periodicity, mirror boundary conditions) can be reset. Further, this
+        also recalculates the appropriate cell size based on the particle
+        kernel radius, `h`. This should be called explicitly when desired but
+        usually this is done when the particles are moved or the `h` is
+        changed.
+
+        The integrator should explicitly call this when needed in the
+        `one_timestep` method.
+        """
+        self.nnps.update_domain()
+
 
 ###############################################################################
 # `EulerIntegrator` class
@@ -244,6 +262,7 @@ class EulerIntegrator(Integrator):
     def one_timestep(self, t, dt):
         self.compute_accelerations()
         self.stage1()
+        self.update_domain()
         self.do_post_stage(dt, 1)
 
 
@@ -269,6 +288,7 @@ def one_timestep(self, t, dt):
 
         # Predict
         self.stage1()
+        self.update_domain()
 
         # Call any post-stage functions.
         self.do_post_stage(0.5*dt, 1)
@@ -277,6 +297,7 @@ def one_timestep(self, t, dt):
 
         # Correct
         self.stage2()
+        self.update_domain()
 
         # Call any post-stage functions.
         self.do_post_stage(dt, 2)
@@ -324,6 +345,7 @@ def one_timestep(self, t, dt):
 
         # Predict
         self.stage1()
+        self.update_domain()
 
         # Call any post-stage functions.
         self.do_post_stage(0.5*dt, 1)
@@ -332,6 +354,7 @@ def one_timestep(self, t, dt):
 
         # Correct
         self.stage2()
+        self.update_domain()
 
         # Call any post-stage functions.
         self.do_post_stage(dt, 2)
@@ -361,16 +384,19 @@ def one_timestep(self, t, dt):
         # stage 1
         self.compute_accelerations()
         self.stage1()
+        self.update_domain()
         self.do_post_stage(1./3*dt, 1)
 
         # stage 2
         self.compute_accelerations()
         self.stage2()
+        self.update_domain()
         self.do_post_stage(2./3*dt, 2)
 
         # stage 3 and end
         self.compute_accelerations()
         self.stage3()
+        self.update_domain()
         self.do_post_stage(dt, 3)
 
 
@@ -382,10 +408,12 @@ class LeapFrogIntegrator(PECIntegrator):
     def one_timestep(self, t, dt):
 
         self.stage1()
+        self.update_domain()
         self.do_post_stage(0.5*dt, 1)
 
         self.compute_accelerations()
         self.stage2()
+        self.update_domain()
         self.do_post_stage(dt, 2)
 
 
@@ -405,20 +433,25 @@ class PEFRLIntegrator(Integrator):
     def one_timestep(self, t, dt):
 
         self.stage1()
+        self.update_domain()
         self.do_post_stage(0.1786178958448091*dt, 1)
 
         self.compute_accelerations()
         self.stage2()
+        self.update_domain()
         self.do_post_stage(0.1123533131749906*dt, 2)
 
         self.compute_accelerations()
         self.stage3()
+        self.update_domain()
         self.do_post_stage(0.8876466868250094*dt, 3)
 
         self.compute_accelerations()
         self.stage4()
+        self.update_domain()
         self.do_post_stage(0.8213821041551909*dt, 4)
 
         self.compute_accelerations()
         self.stage5()
+        self.update_domain()
         self.do_post_stage(dt, 5)

pysph/sph/integrator_cython.mako

@@ -48,6 +48,9 @@ cdef class Integrator:
     cpdef compute_accelerations(self, int index=0, update_nnps=True):
         self.integrator.compute_accelerations(index, update_nnps)
 
+    cpdef update_domain(self):
+        self.integrator.update_domain()
+
     cpdef do_post_stage(self, double stage_dt, int stage):
         """This is called after every stage of the integrator.
 

pysph/sph/integrator_opencl_helper.py

@@ -82,6 +82,9 @@ def set_post_stage_callback(self, callback):
     def compute_accelerations(self, index=0, update_nnps=True):
         self.integrator.compute_accelerations(index, update_nnps)
 
+    def update_domain(self):
+        self.integrator.update_domain()
+
     def do_post_stage(self, stage_dt, stage):
         """This is called after every stage of the integrator.
 

pysph/sph/tests/test_multi_group_integrator.py

@@ -42,6 +42,8 @@ def one_timestep(self, t, dt):
         self.do_post_stage(dt, 1)
         self.compute_accelerations(1, update_nnps=False)
         self.stage2()
+        self.update_domain()
+        self.do_post_stage(dt, 2)
 
 
 class TestMultiGroupIntegrator(unittest.TestCase):
@@ -75,6 +77,7 @@ def _make_integrator(self):
         nnps.update()
         for ae in a_evals:
             ae.set_nnps(nnps)
+        integrator.set_nnps(nnps)
         return integrator
 
     def test_different_accels_per_integrator(self):

pysph/sph/wc/gtvf.py

@@ -52,16 +52,19 @@ def get_particle_array_gtvf(constants=None, **props):
 class GTVFIntegrator(Integrator):
     def one_timestep(self, t, dt):
         self.stage1()
+        self.do_post_stage(dt, 1)
 
         self.compute_accelerations(0, update_nnps=False)
 
         self.stage2()
-
+        # We update domain here alone as positions only change here.
+        self.update_domain()
         self.do_post_stage(dt, 2)
 
         self.compute_accelerations(1)
 
         self.stage3()
+        self.do_post_stage(dt, 3)
 
 
 class GTVFStep(IntegratorStep):

pysph/sph/wc/pcisph.py

@@ -40,8 +40,9 @@ def one_timestep(self, t, dt):
         self.compute_accelerations(0)
 
         self.stage1()
+        self.update_domain()
 
-        self.do_post_stage(0.5 * dt, 1)
+        self.do_post_stage(dt, 1)
 
     def initial_acceleration(self, t, dt):
         pass