Merge pull request #201 from psnbaba/inviscid-wall

added inviscid wall

pysph/examples/flow_past_cylinder_2d.py

@@ -13,7 +13,6 @@
 
 from pysph.base.kernels import QuinticSpline
 from pysph.sph.equation import Equation
-from pysph.base.nnps import DomainManager
 from pysph.base.utils import get_particle_array
 from pysph.solver.application import Application
 from pysph.sph.scheme import SchemeChooser
@@ -40,6 +39,21 @@
 p0 = rho * c0 * c0
 
 
+class ExtrapolateUhat(Equation):
+    def initialize(self, d_idx, d_uhat, d_wij):
+        d_uhat[d_idx] = 0.0
+        d_wij[d_idx] = 0.0
+
+    def loop(self, d_idx, s_idx, d_uhat, s_uhat, d_wij, s_rho,
+             d_au, d_av, d_aw, WIJ, XIJ):
+        d_uhat[d_idx] += s_uhat[s_idx]*WIJ
+        d_wij[d_idx] += WIJ
+
+    def post_loop(self, d_idx, d_wij, d_uhat, d_rho):
+        if d_wij[d_idx] > 1e-14:
+            d_uhat[d_idx] /= d_wij[d_idx]
+
+
 class ResetInletVelocity(Equation):
     def __init__(self, dest, sources, U, V, W):
         self.U = U
@@ -55,16 +69,6 @@ def loop(self, d_idx, d_u, d_v, d_w, d_x, d_y, d_z, d_xn, d_yn, d_zn):
 
 
 class WindTunnel(Application):
-    def create_domain(self):
-        dx = self.dx
-        i_ghost = n_inlet*dx
-        o_ghost = n_outlet*dx
-        domain = DomainManager(
-            xmin=-i_ghost, xmax=l_tunnel+o_ghost, ymin=-w_tunnel,
-            ymax=w_tunnel, periodic_in_y=True
-        )
-        return domain
-
     def add_user_options(self, group):
         group.add_argument(
             "--re", action="store", type=float, dest="re", default=200,
@@ -110,6 +114,39 @@ def _create_fluid(self):
         )
         return fluid
 
+    def _create_wall(self):
+        dx = self.dx
+        h0 = self.hdx * self.dx
+        x0, y0 = np.mgrid[
+            dx/2: l_tunnel+n_inlet*dx+n_outlet*dx: dx, dx/2: n_wall*dx: dx]
+        x0 -= n_inlet*dx
+        y0 -= n_wall*dx+w_tunnel
+        x0 = np.ravel(x0)
+        y0 = np.ravel(y0)
+
+        x1 = np.copy(x0)
+        y1 = np.copy(y0)
+        y1 += n_wall*dx+2*w_tunnel
+        x1 = np.ravel(x1)
+        y1 = np.ravel(y1)
+        x0 = np.concatenate((x0, x1))
+        y0 = np.concatenate((y0, y1))
+        volume = dx*dx
+        wall = get_particle_array(
+            name='wall', x=x0, y=y0, m=volume*rho, rho=rho, h=h0, V=1.0/volume)
+        return wall
+
+    def _set_wall_normal(self, pa):
+        props = ['xn', 'yn', 'zn']
+        for p in props:
+            pa.add_property(p)
+
+        y = pa.y
+        cond = y > 0.0
+        pa.yn[cond] = 1.0
+        cond = y < 0.0
+        pa.yn[cond] = -1.0
+
     def _create_solid(self):
         dx = self.dx
         h0 = self.hdx * self.dx
@@ -159,11 +196,14 @@ def create_particles(self):
         solid = self._create_solid()
         outlet = self._create_outlet()
         inlet = self._create_inlet()
+        wall = self._create_wall()
         G.remove_overlap_particles(fluid, solid, dx_solid=self.dx, dim=2)
 
         particles = [
-            fluid, inlet, outlet, solid
+            fluid, inlet, outlet, solid, wall
         ]
+
+        self._set_wall_normal(wall)
         self.scheme.setup_properties(particles)
 
         return particles
@@ -175,7 +215,8 @@ def create_scheme(self):
 
         edac = EDACScheme(
             ['fluid'], ['solid'], dim=2, rho0=rho, c0=c0, h=h0,
-            pb=None, nu=nu, alpha=0.2, inlet_outlet_manager=self.iom
+            pb=None, nu=nu, alpha=0.2, inlet_outlet_manager=self.iom,
+            inviscid_solids=['wall']
         )
 
         s = SchemeChooser(default='edac', edac=edac)
@@ -197,7 +238,8 @@ def configure_scheme(self):
     def _create_inlet_outlet_manager(self):
         inleteqns = [
             ResetInletVelocity('inlet', [], U=umax, V=0.0, W=0.0),
-            SolidWallPressureBC('inlet', ['fluid'])
+            SolidWallPressureBC('inlet', ['fluid']),
+            ExtrapolateUhat('inlet', ['fluid'])
         ]
 
         inlet_info = InletInfo(

pysph/sph/wc/edac.py

@@ -145,11 +145,10 @@ def __init__(self, dest, sources, gx=0.0, gy=0.0, gz=0.0):
 
         super(SolidWallPressureBC, self).__init__(dest, sources)
 
-    def initialize(self, d_idx, d_p, d_wij):
+    def initialize(self, d_idx, d_p):
         d_p[d_idx] = 0.0
-        d_wij[d_idx] = 0.0
 
-    def loop(self, d_idx, s_idx, d_p, s_p, d_wij, s_rho,
+    def loop(self, d_idx, s_idx, d_p, s_p, s_rho,
              d_au, d_av, d_aw, WIJ, XIJ):
 
         # numerator of Eq. (27) ax, ay and az are the prescribed wall
@@ -161,9 +160,6 @@ def loop(self, d_idx, s_idx, d_p, s_p, d_wij, s_rho,
 
         d_p[d_idx] += s_p[s_idx]*WIJ + s_rho[s_idx]*gdotxij*WIJ
 
-        # denominator of Eq. (27)
-        d_wij[d_idx] += WIJ
-
     def post_loop(self, d_idx, d_wij, d_p):
         # extrapolated pressure at the ghost particle
         if d_wij[d_idx] > 1e-14:
@@ -178,6 +174,15 @@ def post_loop(self, d_idx, d_p):
             d_p[d_idx] = 0.0
 
 
+class SourceNumberDensity(Equation):
+    r"""Evaluates the number density due to the source particles"""
+    def initialize(self, d_idx, d_wij):
+        d_wij[d_idx] = 0.0
+
+    def loop(self, d_idx, d_wij, WIJ):
+        d_wij[d_idx] += WIJ
+
+
 class SetWallVelocity(Equation):
     r"""Extrapolating the fluid velocity on to the wall Eq. (22) in REF1:
 
@@ -195,13 +200,13 @@ class SetWallVelocity(Equation):
     *filtered* velocity variables :math:`uf, vf, wf`.
 
     """
-    def initialize(self, d_idx, d_uf, d_vf, d_wf, d_wij):
+    def initialize(self, d_idx, d_uf, d_vf, d_wf):
         d_uf[d_idx] = 0.0
         d_vf[d_idx] = 0.0
         d_wf[d_idx] = 0.0
 
     def loop(self, d_idx, s_idx, d_uf, d_vf, d_wf,
-             s_u, s_v, s_w, d_wij, WIJ):
+             s_u, s_v, s_w, WIJ):
         # sum in Eq. (22)
         # this will be normalized in post loop
         d_uf[d_idx] += s_u[s_idx] * WIJ
@@ -228,6 +233,71 @@ def post_loop(self, d_uf, d_vf, d_wf, d_wij, d_idx,
         d_wg[d_idx] = 2*d_w[d_idx] - d_wf[d_idx]
 
 
+class NoSlipVelocityExtrapolation(Equation):
+    '''No Slip boundary condition on the wall
+
+    The velocity of the fluid is extrapolated over to the wall using
+    shepard extrapolation. The velocity normal to the wall is reflected back
+    to impose no penetration.
+    '''
+    def initialize(self, d_idx, d_u, d_v, d_w):
+        d_u[d_idx] = 0.0
+        d_v[d_idx] = 0.0
+        d_w[d_idx] = 0.0
+
+    def loop(self, d_idx, s_idx, d_u, d_v, d_w, s_u, s_v, s_w, WIJ,
+             XIJ):
+        d_u[d_idx] += s_u[s_idx]*WIJ
+        d_v[d_idx] += s_v[s_idx]*WIJ
+        d_w[d_idx] += s_w[s_idx]*WIJ
+
+    def post_loop(self, d_idx, d_wij, d_u, d_v, d_w, d_xn, d_yn, d_zn):
+        if d_wij[d_idx] > 1e-14:
+            d_u[d_idx] /= d_wij[d_idx]
+            d_v[d_idx] /= d_wij[d_idx]
+            d_w[d_idx] /= d_wij[d_idx]
+
+        projection = d_u[d_idx]*d_xn[d_idx] +\
+            d_v[d_idx]*d_yn[d_idx] + d_w[d_idx]*d_zn[d_idx]
+
+        d_u[d_idx] = d_u[d_idx] - 2 * projection * d_xn[d_idx]
+        d_v[d_idx] = d_v[d_idx] - 2 * projection * d_yn[d_idx]
+        d_w[d_idx] = d_w[d_idx] - 2 * projection * d_zn[d_idx]
+
+
+class NoSlipAdvVelocityExtrapolation(Equation):
+    '''No Slip boundary condition on the wall
+
+    The advection velocity of the fluid is extrapolated over to the wall
+    using shepard extrapolation. The advection velocity normal to the wall
+    is reflected back to impose no penetration.
+    '''
+    def initialize(self, d_idx, d_uhat, d_vhat, d_what):
+        d_uhat[d_idx] = 0.0
+        d_vhat[d_idx] = 0.0
+        d_what[d_idx] = 0.0
+
+    def loop(self, d_idx, s_idx, d_uhat, d_vhat, d_what, s_uhat, s_vhat,
+             s_what, WIJ, XIJ):
+        d_uhat[d_idx] += s_uhat[s_idx]*WIJ
+        d_vhat[d_idx] += s_vhat[s_idx]*WIJ
+        d_what[d_idx] += s_what[s_idx]*WIJ
+
+    def post_loop(self, d_idx, d_wij, d_uhat, d_vhat, d_what, d_xn, d_yn,
+                  d_zn):
+        if d_wij[d_idx] > 1e-14:
+            d_uhat[d_idx] /= d_wij[d_idx]
+            d_vhat[d_idx] /= d_wij[d_idx]
+            d_what[d_idx] /= d_wij[d_idx]
+
+        projection = d_uhat[d_idx]*d_xn[d_idx] +\
+            d_vhat[d_idx]*d_yn[d_idx] + d_what[d_idx]*d_zn[d_idx]
+
+        d_uhat[d_idx] = d_uhat[d_idx] - 2 * projection * d_xn[d_idx]
+        d_vhat[d_idx] = d_vhat[d_idx] - 2 * projection * d_yn[d_idx]
+        d_what[d_idx] = d_what[d_idx] - 2 * projection * d_zn[d_idx]
+
+
 class MomentumEquation(Equation):
     r"""Momentum equation (gradient of pressure) based on the number
     density formulation of Hu and Adams JCP 213 (2006), 844-861.
@@ -474,7 +544,7 @@ class EDACScheme(Scheme):
     def __init__(self, fluids, solids, dim, c0, nu, rho0, pb=0.0,
                  gx=0.0, gy=0.0, gz=0.0, tdamp=0.0, eps=0.0, h=0.0,
                  edac_alpha=0.5, alpha=0.0, bql=True, clamp_p=False,
-                 inlet_outlet_manager=None):
+                 inlet_outlet_manager=None, inviscid_solids=None):
         """The EDAC scheme.
 
         Parameters
@@ -512,6 +582,10 @@ def __init__(self, fluids, solids, dim, c0, nu, rho0, pb=0.0,
         clamp_p : bool
             Clamp the boundary pressure to positive values.  This is only used
             for external flows.
+        inlet_outlet_manager : InletOutletManager Instance
+            Pass the manager if inlet outlet boundaries are present
+        inviscid_solids : list
+            list of inviscid solid array names
         """
         self.c0 = c0
         self.nu = nu
@@ -532,6 +606,8 @@ def __init__(self, fluids, solids, dim, c0, nu, rho0, pb=0.0,
         self.alpha = alpha
         self.h = h
         self.inlet_outlet_manager = inlet_outlet_manager
+        self.inviscid_solids = [] if inviscid_solids is None else\
+            inviscid_solids
         self.attributes_changed()
 
     # Public protocol ###################################################
@@ -673,11 +749,14 @@ def setup_properties(self, particles, clean=True):
             if iom is not None:
                 iom.add_io_properties(pa, self)
 
-        TVF_SOLID_PROPS = ['V', 'wij', 'ax', 'ay', 'az',
-                           'uf', 'vf', 'wf', 'ug', 'vg', 'wg']
+        TVF_SOLID_PROPS = ['V', 'wij', 'ax', 'ay', 'az', 'uf', 'vf', 'wf',
+                           'ug', 'vg', 'wg']
+        if self.inviscid_solids:
+            TVF_SOLID_PROPS += ['xn', 'yn', 'zn', 'uhat',
+                                'vhat', 'what']
         extra_props = TVF_SOLID_PROPS if self.use_tvf else EDAC_SOLID_PROPS
         all_solid_props = DEFAULT_PROPS.union(extra_props)
-        for solid in self.solids:
+        for solid in (self.solids+self.inviscid_solids):
             pa = particle_arrays[solid]
             self._ensure_properties(pa, all_solid_props, clean)
             pa.set_output_arrays(['x', 'y', 'z', 'u', 'v', 'w', 'rho', 'p',
@@ -704,9 +783,10 @@ def _get_internal_flow_equations(self):
 
         iom = self.inlet_outlet_manager
         fluids_with_io = self.fluids
+        all_solids = self.solids + self.inviscid_solids
         if iom is not None:
             fluids_with_io = self.fluids + iom.get_io_names()
-        all = fluids_with_io + self.solids
+        all = fluids_with_io + all_solids
 
         equations = []
         # inlet-outlet
@@ -717,7 +797,7 @@ def _get_internal_flow_equations(self):
 
         group1 = []
         avg_p_group = []
-        has_solids = len(self.solids) > 0
+        has_solids = len(all_solids) > 0
         for fluid in fluids_with_io:
             group1.append(SummationDensity(dest=fluid, sources=all))
             if self.bql:
@@ -729,11 +809,23 @@ def _get_internal_flow_equations(self):
 
         for solid in self.solids:
             group1.extend([
+                SourceNumberDensity(dest=solid, sources=fluids_with_io),
                 VolumeSummation(dest=solid, sources=all),
                 SolidWallPressureBC(dest=solid, sources=fluids_with_io,
                                     gx=self.gx, gy=self.gy, gz=self.gz),
                 SetWallVelocity(dest=solid, sources=fluids_with_io),
             ])
+        for solid in self.inviscid_solids:
+            group1.extend([
+                SourceNumberDensity(dest=solid, sources=fluids_with_io),
+                NoSlipVelocityExtrapolation(
+                    dest=solid, sources=fluids_with_io),
+                NoSlipAdvVelocityExtrapolation(
+                    dest=solid, sources=fluids_with_io),
+                VolumeSummation(dest=solid, sources=all),
+                SolidWallPressureBC(dest=solid, sources=fluids_with_io,
+                                    gx=self.gx, gy=self.gy, gz=self.gz)
+                ])
 
         equations.append(Group(equations=group1, real=False))
 
@@ -780,6 +872,7 @@ def _get_internal_flow_equations(self):
             ])
         equations.append(Group(equations=group2))
 
+        print(equations)
         return equations
 
     def _get_external_flow_equations(self):
@@ -790,7 +883,8 @@ def _get_external_flow_equations(self):
             MomentumEquationArtificialViscosity,
             MomentumEquationViscosity
         )
-        all = self.fluids + self.solids
+        all_solids = self.solids + self.inviscid_solids
+        all = self.fluids + all_solids
 
         edac_nu = self._get_edac_nu()
         equations = []
@@ -799,6 +893,7 @@ def _get_external_flow_equations(self):
             group1.append(SummationDensity(dest=fluid, sources=all))
         for solid in self.solids:
             group1.extend([
+                SourceNumberDensity(dest=solid, sources=self.fluids),
                 VolumeSummation(dest=solid, sources=all),
                 SolidWallPressureBC(dest=solid, sources=self.fluids,
                                     gx=self.gx, gy=self.gy, gz=self.gz),
@@ -809,6 +904,16 @@ def _get_external_flow_equations(self):
                     ClampWallPressure(dest=solid, sources=None)
                 )
 
+        for solid in self.inviscid_solids:
+            group1.extend([
+                SourceNumberDensity(dest=solid, sources=self.fluids),
+                NoSlipVelocityExtrapolation(
+                    dest=solid, sources=self.fluids),
+                VolumeSummation(dest=solid, sources=all),
+                SolidWallPressureBC(dest=solid, sources=self.fluids,
+                                    gx=self.gx, gy=self.gy, gz=self.gz)
+                ])
+
         equations.append(Group(equations=group1, real=False))
 
         group2 = []