Merge pull request #232 from abhinavmuta/sisph

Simple iterative incompressible SPH

pysph/examples/dam_break_2d.py

@@ -11,7 +11,7 @@
 import numpy as np
 
 from pysph.base.kernels import WendlandQuintic, QuinticSpline
-from pysph.examples._db_geometry import DamBreak2DGeometry
+from pysph.base.utils import get_particle_array
 from pysph.solver.application import Application
 from pysph.sph.scheme import WCSPHScheme, SchemeChooser, AdamiHuAdamsScheme
 from pysph.sph.wc.edac import EDACScheme
@@ -22,6 +22,8 @@
                                             MixedKernelCorrectionPreStep)
 from pysph.sph.wc.crksph import CRKSPHPreStep, CRKSPH
 from pysph.sph.wc.gtvf import GTVFScheme
+from pysph.sph.isph.sisph import SISPHScheme
+from pysph.tools.geometry import get_2d_tank, get_2d_block
 
 
 fluid_column_height = 2.0
@@ -33,14 +35,16 @@
 dx = 0.03
 g = 9.81
 ro = 1000.0
-co = 10.0 * np.sqrt(2 * 9.81 * fluid_column_height)
+vref = np.sqrt(2 * 9.81 * fluid_column_height)
+co = 10.0 * vref
 gamma = 7.0
 alpha = 0.1
 beta = 0.0
 B = co * co * ro / gamma
 p0 = 1000.0
 hdx = 1.3
 h = hdx * dx
+m = dx**2 * ro
 
 
 class DamBreak2D(Application):
@@ -124,6 +128,14 @@ def configure_scheme(self):
             kw.update(dict(kernel=kernel, dt=dt))
             scheme.configure(pref=B*gamma, h0=self.h)
             print("dt = %f" % dt)
+        elif self.options.scheme == 'sisph':
+            dt = 0.125*self.h/vref
+            kernel = QuinticSpline(dim=2)
+            print("SISPH dt = %f" % dt)
+            kw.update(dict(kernel=kernel))
+            self.scheme.configure_solver(
+                dt=dt, tf=tf, adaptive_timestep=False, pfreq=10,
+            )
 
         self.scheme.configure_solver(**kw)
 
@@ -149,13 +161,18 @@ def create_scheme(self):
             fluids=['fluid'], solids=['boundary'], dim=2, nu=nu,
             rho0=ro, gy=-g, h0=None, c0=co, pref=None
         )
+        sisph = SISPHScheme(
+            fluids=['fluid'], solids=['boundary'], dim=2, nu=nu,
+            c0=co, rho0=ro, alpha=0.05, gy=-g, pref=ro*co**2,
+            internal_flow=False, hg_correction=True, gtvf=True, symmetric=True
+        )
         s = SchemeChooser(default='wcsph', wcsph=wcsph, aha=aha, edac=edac,
-                          iisph=iisph, gtvf=gtvf)
+                          iisph=iisph, gtvf=gtvf, sisph=sisph)
         return s
 
     def create_equations(self):
         eqns = self.scheme.get_equations()
-        if self.options.scheme == 'iisph':
+        if self.options.scheme == 'iisph' or self.options.scheme == 'sisph':
             return eqns
         if self.options.scheme == 'gtvf':
             return eqns
@@ -198,17 +215,19 @@ def create_particles(self):
             nboundary_layers = 4
             nfluid_offset = 1
             wall_hex_pack = False
-        geom = DamBreak2DGeometry(
-            container_width=container_width, container_height=container_height,
-            fluid_column_height=fluid_column_height,
-            fluid_column_width=fluid_column_width, dx=self.dx, dy=self.dx,
-            nboundary_layers=1, ro=ro, co=co,
-            with_obstacle=False, wall_hex_pack=wall_hex_pack,
-            beta=1.0, nfluid_offset=1, hdx=self.hdx)
-        fluid, boundary = geom.create_particles(
-            nboundary_layers=nboundary_layers, hdx=self.hdx,
-            nfluid_offset=nfluid_offset,
-        )
+        xt, yt = get_2d_tank(dx=self.dx, length=container_width,
+                             height=container_height, base_center=[2, 0],
+                             num_layers=nboundary_layers)
+        xf, yf = get_2d_block(dx=self.dx, length=fluid_column_width,
+                              height=fluid_column_height, center=[0.5, 1])
+
+        xf += self.dx
+        yf += self.dx
+
+        fluid = get_particle_array(name='fluid', x=xf, y=yf, h=h, m=m, rho=ro)
+        boundary = get_particle_array(name='boundary', x=xt, y=yt, h=h, m=m,
+                                      rho=ro)
+
         self.scheme.setup_properties([fluid, boundary])
         if self.options.scheme == 'iisph':
             # the default position tends to cause the particles to be pushed
@@ -267,6 +286,12 @@ def post_process(self, info_fname):
         plt.savefig(os.path.join(self.output_dir, 'x_vs_t.png'), dpi=300)
         plt.close()
 
+    def customize_output(self):
+        self._mayavi_config('''
+        b = particle_arrays['fluid']
+        b.scalar = 'vmag'
+        ''')
+
 
 if __name__ == '__main__':
     app = DamBreak2D()

pysph/examples/taylor_green.py

@@ -1,9 +1,9 @@
 """Taylor Green vortex flow (5 minutes).
 """
 
+import os
 import numpy as np
 from numpy import pi, sin, cos, exp
-import os
 
 from pysph.base.nnps import DomainManager
 from pysph.base.utils import get_particle_array
@@ -23,6 +23,7 @@
 from pysph.sph.wc.gtvf import GTVFScheme
 from pysph.sph.wc.pcisph import PCISPHScheme
 from pysph.sph.wc.shift import ShiftPositions
+from pysph.sph.isph.sisph import SISPHScheme
 
 
 # domain and constants
@@ -152,6 +153,8 @@ def consume_user_options(self):
         h0 = self.hdx * self.dx
         if self.options.scheme == 'iisph' or self.options.scheme == 'pcisph':
             dt_cfl = 0.25 * h0 / U
+        if self.options.scheme == 'sisph':
+            dt_cfl = 0.25 * h0 / U
         else:
             dt_cfl = 0.25 * h0 / (c0 + U)
         dt_viscous = 0.125 * h0**2 / nu
@@ -179,6 +182,9 @@ def configure_scheme(self):
             scheme.configure(h0=h0, nu=self.nu)
         elif self.options.scheme == 'gtvf':
             scheme.configure(pref=p0, nu=self.nu, h0=h0)
+        elif self.options.scheme == 'sisph':
+            pfreq = 10
+            scheme.configure(nu=self.nu)
         scheme.configure_solver(kernel=kernel, tf=self.tf, dt=self.dt,
                                 pfreq=pfreq)
 
@@ -212,9 +218,14 @@ def create_scheme(self):
         pcisph = PCISPHScheme(
             fluids=['fluid'], dim=2, rho0=rho0, nu=None
         )
+        sisph = SISPHScheme(
+            fluids=['fluid'], solids=[], dim=2, nu=None, rho0=rho0,
+            c0=c0, alpha=0.0, has_ghosts=True, pref=p0,
+            rho_cutoff=0.2, internal_flow=True, gtvf=True
+        )
         s = SchemeChooser(
             default='tvf', wcsph=wcsph, tvf=tvf, edac=edac, iisph=iisph,
-            crksph=crksph, gtvf=gtvf, pcisph=pcisph
+            crksph=crksph, gtvf=gtvf, pcisph=pcisph, sisph=sisph
         )
         return s
 
@@ -290,6 +301,10 @@ def create_particles(self):
             fluid.get_number_of_particles(), self.dt))
 
         # volume is set as dx^2
+        if self.options.scheme == 'sisph':
+            nfp = fluid.get_number_of_particles()
+            fluid.gid[:] = np.arange(nfp)
+            fluid.add_output_arrays(['gid'])
         if self.options.scheme == 'tvf':
             fluid.V[:] = 1. / self.volume
         if self.options.scheme == 'iisph':
@@ -483,6 +498,12 @@ def post_process(self, info_fname):
         fig = os.path.join(self.output_dir, "p_l1_error.png")
         plt.savefig(fig, dpi=300)
 
+    def customize_output(self):
+        self._mayavi_config('''
+        b = particle_arrays['fluid']
+        b.scalar = 'vmag'
+        ''')
+
 
 if __name__ == '__main__':
     app = TaylorGreen()