Merged in prabhuramachandran/pysph/fix_iisph (pull request #141)

Fix some issues with IISPH implementation.
pypr · Feb 1, 2015 · 8327c7e · 8327c7e
2 parents 159e0fc + 120040c
commit 8327c7e
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 22 deletions.
diff --git a/examples/iisph/dam_break.py b/examples/iisph/dam_break.py
@@ -72,6 +72,7 @@
 
 """
 
+import numpy as np
 import sys
 import os
 # Need this to import db_geometry.
@@ -101,14 +102,15 @@
 container_width  = 4.0
 nboundary_layers=2
 
-dt = 1e-2
+dt = 1e-3
 tf = 2.5
 
-hdx = 1.0
+hdx = 1.3
 dx = dy = 0.02
 ro = 1000.0
 nu = 8.9e-4
 beta = 0.0
+co = 10.0*np.sqrt(9.81*fluid_column_height)
 
 geom = DamBreak2DGeometry(
     container_width=container_width, container_height=container_height,
@@ -121,8 +123,8 @@
 
 def create_particles(**kw):
     fluid, boundary = geom.create_particles(**kw)
-    boundary.x -= 0.025
-    boundary.y -= 0.025
+    boundary.x -= 0.1
+    boundary.y -= 0.1
     return [fluid, boundary]
 
 # Create the application.
@@ -177,7 +179,7 @@ def create_particles(**kw):
                 dest='fluid', sources=['boundary'], nu=nu, rho0=ro,
             ),
             ComputeDII(dest='fluid', sources=['fluid']),
-            ComputeDIIBoundary(dest='fluid', sources=['boundary'], rho0=ro),
+            #ComputeDIIBoundary(dest='fluid', sources=['boundary'], rho0=ro),
         ]
     ),
 
@@ -186,7 +188,7 @@ def create_particles(**kw):
             ComputeRhoAdvection(dest='fluid', sources=['fluid']),
             ComputeRhoBoundary(dest='fluid', sources=['boundary'], rho0=ro),
             ComputeAII(dest='fluid', sources=['fluid']),
-            ComputeAIIBoundary(dest='fluid', sources=['boundary'], rho0=ro),
+            #ComputeAIIBoundary(dest='fluid', sources=['boundary'], rho0=ro),
         ]
     ),
 
@@ -203,14 +205,14 @@ def create_particles(**kw):
                         dest='fluid', sources=['fluid'], rho0=ro,
                         tolerance=1e-3, debug=False
                     ),
-                    PressureSolveBoundary(
-                        dest='fluid', sources=['boundary'], rho0=ro,
-                    ),
+                    #PressureSolveBoundary(
+                    #    dest='fluid', sources=['boundary'], rho0=ro,
+                    #),
                   ]
             ),
         ],
         iterate=True,
-        max_iterations=20,
+        max_iterations=30,
         min_iterations=2
     ),
 

diff --git a/pysph/sph/iisph.py b/pysph/sph/iisph.py
@@ -101,11 +101,11 @@ def __init__(self, dest, sources=None, nu=1.0):
 
     def loop(self, d_idx, d_au, d_av, d_aw, s_idx, s_m, EPS,
              VIJ, XIJ, RHOIJ1, R2IJ, DWIJ):
-        vijdotxij = min(VIJ[0]*XIJ[0] + VIJ[1]*XIJ[1] + VIJ[2]*XIJ[2], 0.0)
-        fac = self.nu*s_m[s_idx]*RHOIJ1*vijdotxij/(R2IJ + EPS)
-        d_au[d_idx] += fac*DWIJ[0]
-        d_av[d_idx] += fac*DWIJ[1]
-        d_aw[d_idx] += fac*DWIJ[2]
+        dwijdotxij = DWIJ[0]*XIJ[0] + DWIJ[1]*XIJ[1] + DWIJ[2]*XIJ[2]
+        fac = 2.0*self.nu*s_m[s_idx]*RHOIJ1*dwijdotxij/(R2IJ + EPS)
+        d_au[d_idx] += fac*VIJ[0]
+        d_av[d_idx] += fac*VIJ[1]
+        d_aw[d_idx] += fac*VIJ[2]
 
 class ViscosityAccelerationBoundary(Equation):
     """The acceleration on the fluid due to a boundary.
@@ -118,12 +118,12 @@ def __init__(self, dest, sources=None, rho0=1.0, nu=1.0):
     def loop(self, d_idx, d_au, d_av, d_aw, d_rho, s_idx, s_V, EPS,
              VIJ, XIJ, R2IJ, DWIJ):
         phi_b = self.rho0/(s_V[s_idx]*d_rho[d_idx])
-        vijdotxij = min(VIJ[0]*XIJ[0] + VIJ[1]*XIJ[1] + VIJ[2]*XIJ[2], 0.0)
+        dwijdotxij = DWIJ[0]*XIJ[0] + DWIJ[1]*XIJ[1] + DWIJ[2]*XIJ[2]
 
-        fac = self.nu*phi_b*vijdotxij/(R2IJ + EPS)
-        d_au[d_idx] += fac*DWIJ[0]
-        d_av[d_idx] += fac*DWIJ[1]
-        d_aw[d_idx] += fac*DWIJ[2]
+        fac = 2.0*self.nu*phi_b*dwijdotxij/(R2IJ + EPS)
+        d_au[d_idx] += fac*VIJ[0]
+        d_av[d_idx] += fac*VIJ[1]
+        d_aw[d_idx] += fac*VIJ[2]
 
 
 class ComputeDII(Equation):
@@ -339,8 +339,8 @@ def loop(self, d_idx, d_rho, d_p, d_au, d_av, d_aw,
 
     def post_loop(self, d_idx, d_au, d_av, d_aw,
                   d_uadv, d_vadv, d_wadv, DT_ADAPT):
-        fac = sqrt(d_au[d_idx]*d_au[d_idx] + d_av[d_idx]*d_av[d_idx] +\
-                   d_aw[d_idx]*d_aw[d_idx])
+        fac = d_au[d_idx]*d_au[d_idx] + d_av[d_idx]*d_av[d_idx] +\
+                   d_aw[d_idx]*d_aw[d_idx]
         vmag = sqrt(d_uadv[d_idx]*d_uadv[d_idx] + d_vadv[d_idx]*d_vadv[d_idx] +
                     d_wadv[d_idx]*d_wadv[d_idx])
         DT_ADAPT[0] = max(2.0*vmag, DT_ADAPT[0])