Merged in prabhuramachandran/pysph/inlet_outlet (pull request #145)

Adding support for simple inlets and outlets

examples/trivial_inlet_outlet.py

@@ -0,0 +1,110 @@
+"""A trivial example to demonstrate the inlet and outlet feature in 2D.
+
+We first create three particle arrays, an "inlet", "fluid" and "outlet" in the
+`create_particle` function. Initially there are no fluid and outlet particles.
+A single row of inlet particles are created between (0.0, 0.0) to (0.0, 1.0),
+i.e. along the y-axis with a u velocity = 0.25.
+
+The inlet and outlet are created in the `create_inlet_outlet` function.  This
+function is passed a dictionary of `{array_name:particle_array}`. An inlet
+between (-0.4, 0.0) and (0.0, 1.0) is created by instantiating a
+`SimpleInlet`.  The inlet first makes 4 copies of the inlet particle array
+data and stacks them along the negative x-axis.  The `InletOutletStep` is used
+to step all particles and simply moves the particles.  As particles leave
+the inlet they are converted to fluid particles.
+
+An outlet is also created in the region (0.5, 0.0), (1.0, 1.0) and as fluid
+particles enter the outlet region, they are converted to outlet particles.  As
+outlet particles leave the outlet they are removed from the simulation.
+
+The following figure should make this clear.
+
+               inlet       fluid       outlet
+              ---------    --------    --------
+             | * * * x |  |        |  |        |
+     u       | * * * x |  |        |  |        |
+    --->     | * * * x |  |        |  |        |
+             | * * * x |  |        |  |        |
+              --------     --------    --------
+
+In the figure above, the 'x' are the initial inlet particles.  The '*' are the
+copies of these.  The particles are moving to the right and as they do, new
+fluid particles are added and as the fluid particles flow into the outlet they
+are converted to the outlet particle array and at last as the particles leave
+the outlet they are removed from the simulation.  The `create_particles` and
+`create_inlet_outlet` functions are passed to the `app.setup` method.
+
+This example can be run in parallel.
+
+"""
+
+import numpy as np
+
+from pysph.base.kernels import CubicSpline
+from pysph.base.utils import get_particle_array
+from pysph.solver.application import Application
+from pysph.solver.solver import Solver
+from pysph.sph.integrator import PECIntegrator
+from pysph.sph.simple_inlet_outlet import SimpleInlet, SimpleOutlet
+from pysph.sph.integrator_step import InletOutletStep
+
+def create_particles():
+    # Note that you need to create the inlet and outlet arrays in this method.
+
+    # Initially fluid has no particles -- these are generated by the inlet.
+    fluid = get_particle_array(name='fluid')
+
+    outlet = get_particle_array(name='outlet')
+
+    # Setup the inlet particle array with just the particles we need at the
+    # exit plane which is replicated by the inlet.
+    dx = 0.1
+    y = np.linspace(0, 1, 11)
+    x = np.zeros_like(y)
+    m = np.ones_like(x)*dx
+    h = np.ones_like(x)*dx*1.5
+    rho = np.ones_like(x)
+    # Remember to set u otherwise the inlet particles won't move.
+    u = np.ones_like(x)*0.25
+
+    inlet = get_particle_array(name='inlet', x=x, y=y, m=m, h=h, u=u, rho=rho)
+
+    return [inlet, fluid, outlet]
+
+def create_inlet_outlet(particle_arrays):
+    # particle_arrays is a dict {name: particle_array}
+    fluid_pa = particle_arrays['fluid']
+    inlet_pa = particle_arrays['inlet']
+    outlet_pa = particle_arrays['outlet']
+
+    # Create the inlet and outlets as described in the documentation.
+    inlet = SimpleInlet(
+        inlet_pa, fluid_pa, spacing=0.1, n=5, axis='x', xmin=-0.4, xmax=0.0,
+        ymin=0.0, ymax=1.0
+    )
+    outlet = SimpleOutlet(
+        outlet_pa, fluid_pa, xmin=0.5, xmax=1.0, ymin=0.0, ymax=1.0
+    )
+    return [inlet, outlet]
+
+
+app = Application()
+kernel = CubicSpline(dim=2)
+integrator = PECIntegrator(
+    fluid=InletOutletStep(), inlet=InletOutletStep(), outlet=InletOutletStep()
+)
+
+dt = 1e-2
+tf = 6
+
+solver = Solver(
+    kernel=kernel, dim=2, integrator=integrator, dt=dt, tf=tf,
+    adaptive_timestep=False
+)
+
+app.setup(
+    solver=solver, equations=[], particle_factory=create_particles,
+    inlet_outlet_factory=create_inlet_outlet
+)
+
+app.run()

pysph/parallel/parallel_manager.pyx

@@ -219,6 +219,7 @@ cdef class ParticleArrayExchange:
             zcomm.Comm_Do( sendbuf, recvbuf )
 
     def remove_remote_particles(self):
+        self.num_local = self.pa.get_number_of_particles(real=True)
         cdef int num_local = self.num_local
 
         # resize the particle array
@@ -419,13 +420,13 @@ cdef class ParallelManager:
         """Peform a reduction to compute the globally stable time steps"""
         cdef np.ndarray dt_sendbuf = self.dt_sendbuf
         cdef np.ndarray dt_recvbuf = np.zeros_like(dt_sendbuf)
-        
+
         comm = self.comm
 
         # set the local time step and peform the global reduction
         dt_sendbuf[0] = local_dt
         comm.Allreduce( sendbuf=dt_sendbuf, recvbuf=dt_recvbuf, op=mpi.MIN )
-        
+
         return dt_recvbuf[0]
 
     cpdef compute_cell_size(self):
@@ -1284,15 +1285,15 @@ cdef class ZoltanParallelManagerGeometric(ZoltanParallelManager):
 
     def set_zoltan_rcb_directions(self, value):
         """Option to group the cuts along a given direction
-        
+
         Legal values (refer to the Zoltan User Guide):
-        
-        '0' = don't order cuts; 
-        '1' = xyz 
-        '2' = xzy 
-        '3' = yzx 
-        '4' = yxz 
-        '5' = zxy 
+
+        '0' = don't order cuts;
+        '1' = xyz
+        '2' = xzy
+        '3' = yzx
+        '4' = yxz
+        '5' = zxy
         '6' = zyx
 
         """

pysph/solver/application.py

@@ -77,6 +77,8 @@ def __init__(self, fname=None, domain=None):
         self._setup_optparse()
 
         self.path = None
+        self.particles = []
+        self.inlet_outlet = []
 
     def _setup_optparse(self):
         usage = """
@@ -407,27 +409,38 @@ def _setup_logging(self, filename=None, loglevel=logging.WARNING,
         if stream:
             logger.addHandler(logging.StreamHandler())
 
+    def _create_inlet_outlet(self, inlet_outlet_factory):
+        """Create the inlets and outlets if needed.
+
+        This method requires that the particles be already created.
+
+        The `inlet_outlet_factory` is passed a dictionary of the particle
+        arrays.  The factory should return a list of inlets and outlets.
+        """
+        if inlet_outlet_factory is not None:
+            solver = self._solver
+            particle_arrays = dict([(p.name, p) for p in self.particles])
+            self.inlet_outlet = inlet_outlet_factory(particle_arrays)
+            # Hook up the inlet/outlet's update method to be called after
+            # each stage.
+            for obj in self.inlet_outlet:
+                solver.add_post_step_callback(obj.update)
+
     def _create_particles(self, particle_factory, *args, **kw):
 
         """ Create particles given a callable `particle_factory` and any
         arguments to it.
-
-        This will also automatically distribute the particles among
-        processors if this is a parallel run.  Returns a list of particle
-        arrays that are created.
         """
-
         solver = self._solver
-        num_procs = self.num_procs
         options = self.options
         rank = self.rank
-        comm = self.comm
 
         # particle array info that is used to create dummy particles
         # on non-root processors
         particles_info = {}
 
-        # Only master creates the particles.
+        # Only master actually calls the particle factory, the rest create
+        # dummy particle arrays.
         if rank == 0:
             if options.restart_file is not None:
                 data = load(options.restart_file)
@@ -464,7 +477,16 @@ def _create_particles(self, particle_factory, *args, **kw):
         if rank != 0:
             self.particles = utils.create_dummy_particles(particles_info)
 
+    def _do_initial_load_balancing(self):
+        """ This will automatically distribute the particles among processors
+        if this is a parallel run.
+        """
         # Instantiate the Parallel Manager here and do an initial LB
+        num_procs = self.num_procs
+        options = self.options
+        solver = self._solver
+        comm = self.comm
+
         self.pm = None
         if num_procs > 1:
             options = self.options
@@ -538,8 +560,6 @@ def _create_particles(self, particle_factory, *args, **kw):
         # set the solver's parallel manager
         solver.set_parallel_manager(self.pm)
 
-        return self.particles
-
     ######################################################################
     # Public interface.
     ######################################################################
@@ -557,8 +577,8 @@ def add_option(self, opt):
             for o in opt:
                 self.add_option(o)
 
-    def setup(self, solver, equations, nnps=None, particle_factory=None,
-              *args, **kwargs):
+    def setup(self, solver, equations, nnps=None, inlet_outlet_factory=None,
+              particle_factory=None, *args, **kwargs):
         """Set the application's solver.  This will call the solver's
         `setup` method.
 
@@ -576,17 +596,28 @@ def setup(self, solver, equations, nnps=None, particle_factory=None,
 
         dir -- the output directory
 
-        integration_type -- The integration method
-
-        default_kernel -- the default kernel to use for operations
-
         Parameters
         ----------
+        solver: The solver instance.
+
+        equations: list: a list of Groups/Equations.
+
+        nnps: NNPS instance: optional NNPS instance.
+            If none is given a default NNPS is created.
+
+        inlet_outlet_factory: callable or None
+            The `inlet_outlet_factory` is passed a dictionary of the particle
+            arrays.  The factory should return a list of inlets and outlets.
+
         particle_factory : callable or None
             If supplied, particles will be created for the solver using the
             particle arrays returned by the callable. Else particles for the
             solver need to be set before calling this method
 
+        args: extra arguments passed on to the particle_factory.
+
+        kwargs: extra keyword arguments passed to the particle factory.
+
         """
         start_time = time.time()
         self._solver = solver
@@ -600,6 +631,12 @@ def setup(self, solver, equations, nnps=None, particle_factory=None,
         # Create particles either from scratch or restart
         self._create_particles(particle_factory, *args, **kwargs)
 
+        # This must be done before the initial load balancing
+        # as the inlets will create new particles.
+        self._create_inlet_outlet(inlet_outlet_factory)
+
+        self._do_initial_load_balancing()
+
         # setup the solver using any options
         self._solver.setup_solver(options.__dict__)
 

pysph/sph/acceleration_eval_cython.mako

@@ -256,6 +256,7 @@ cdef class AccelerationEval:
         ## all_eqs is a Group of all equations having this destination.
         #######################################################################
         % for g_idx, group in enumerate(helper.object.mega_groups):
+        % if len(group.data) > 0: # No equations in this group.
         # ---------------------------------------------------------------------
         # Group ${g_idx}.
         % if group.iterate:
@@ -289,5 +290,6 @@ cdef class AccelerationEval:
 
         # Group ${g_idx} done.
         # ---------------------------------------------------------------------
+        % endif # (if len(group.data) > 0)
         % endfor
         self._set_dt_adapt(DT_ADAPT)

pysph/sph/integrator_step.py

@@ -584,7 +584,7 @@ def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, dt=0.0):
         d_v[d_idx] += dtb2 * d_av[d_idx]
         d_w[d_idx] += dtb2 * d_aw[d_idx]
 
-    def stage2(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, 
+    def stage2(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w,
                d_au, d_av, d_aw, dt=0.0):
 
         dtb2 = 0.5 * dt
@@ -594,7 +594,28 @@ def stage2(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w,
         d_y[d_idx] += dt * d_v[d_idx]
         d_z[d_idx] += dt * d_w[d_idx]
 
-        # Eq. (5.53) in [JM05] 
+        # Eq. (5.53) in [JM05]
         d_u[d_idx] += dtb2 * d_au[d_idx]
         d_v[d_idx] += dtb2 * d_av[d_idx]
         d_w[d_idx] += dtb2 * d_aw[d_idx]
+
+###############################################################################
+# `InletOutletStep` class
+###############################################################################
+class InletOutletStep(IntegratorStep):
+    """A trivial integrator for the inlet/outlet particles
+    """
+    def initialize(self):
+        pass
+
+    def stage1(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, dt=0.0):
+        dtb2 = 0.5*dt
+        d_x[d_idx] += dtb2 * d_u[d_idx]
+        d_y[d_idx] += dtb2 * d_v[d_idx]
+        d_z[d_idx] += dtb2 * d_w[d_idx]
+
+    def stage2(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, dt=0.0):
+        dtb2 = 0.5*dt
+        d_x[d_idx] += dtb2 * d_u[d_idx]
+        d_y[d_idx] += dtb2 * d_v[d_idx]
+        d_z[d_idx] += dtb2 * d_w[d_idx]