working on widget need to add controlling fluid box

WCSPH.py

@@ -7,7 +7,8 @@ def __init__(self, particle_system):
         super().__init__(particle_system)
         self.gamma = self.ps.config['gamma']
         self.B = self.ps.config['B']
-        self.surface_tension = 0.01
+        self.surface_tension = ti.field(ti.f32, shape=())
+        self.surface_tension[None] = self.ps.config['surfaceTension']
 
     @ti.func
     def update_density_task(self, p_i, p_j, density: ti.template()):
@@ -93,9 +94,9 @@ def compute_non_pressure_force_task(self, p_i, p_j, acc: ti.template()):
         # One thing to notice is that there should be extra term (x_a-x_b) in (16) check memo in the paper.
         if self.ps.material[p_j] == self.ps.material_fluid:
             r_vec = self.ps.position[p_i] - self.ps.position[p_j]
-            #if r_vec.norm() > self.ps.particle_diameter:
-            acc -= self.surface_tension / self.ps.mass[p_i] * self.ps.mass[p_j] * r_vec * self.cubic_spline_kernel(
-                r_vec.norm())
+            # if r_vec.norm() > self.ps.particle_diameter:
+            acc -= self.surface_tension[None] / self.ps.mass[p_i] * self.ps.mass[p_j] * r_vec * \
+                   self.cubic_spline_kernel(r_vec.norm())
             """
             else:
                 acc -= self.surface_tension / self.ps.mass[p_i] * self.ps.mass[p_j] * r_vec * self.cubic_spline_kernel(
@@ -105,7 +106,8 @@ def compute_non_pressure_force_task(self, p_i, p_j, acc: ti.template()):
         # Viscosity Force
         # Versatile Rigid-Fluid Coupling for Incompressible SPH  equation (11)~(14)
         if self.ps.material[p_j] == self.ps.material_fluid:
-            nu = 2 * self.viscosity[None] * self.ps.support_length * self.c_s / (self.ps.density[p_i] + self.ps.density[p_j])
+            nu = 2 * self.viscosity[None] * self.ps.support_length * self.c_s / (
+                    self.ps.density[p_i] + self.ps.density[p_j])
             v_ij = self.ps.velocity[p_i] - self.ps.velocity[p_j]
             x_ij = self.ps.position[p_i] - self.ps.position[p_j]
             pi = -nu * ti.min(v_ij.dot(x_ij), 0.0) / (x_ij.dot(x_ij) + 0.01 * self.ps.support_length ** 2)  # eq (11)

data/scenes/dragon_bath.json

@@ -17,6 +17,7 @@
 		"dt": 4e-4,
 		"collisionFactor": 0.5,
 		"viscosity": 0.01,
+		"surfaceTension": 0.01,
 		"c_s": 88.5
 	},
 	"RigidBodies": [

run_simulation.py

@@ -71,6 +71,22 @@
     solver.dt[None] = gui.slider_float('[10^-3]', solver.dt[None] * 1000, 0.2, 0.8) * 0.001
     gui.text('Viscosity')
     solver.viscosity[None] = gui.slider_float('', solver.viscosity[None], 0.001, 0.5)
+    gui.text('Surface Tension')
+    solver.surface_tension[None] = gui.slider_float('[N/m]', solver.surface_tension[None], 0.001, 5)
+    if solver.viscosity[None] > 0.23 or solver.surface_tension[None] > 2.0:
+        # Viscosity with over 0.23 cause numerical instability when time step is larger than 0.0005 typically.
+        # Surface tension with over 2.0 cause numerical instability when time step is larger than 0.0005 typically.
+        solver.dt[None] = ti.min(solver.dt[None], 0.0005)
+    if solver.viscosity[None] > 0.23 and solver.surface_tension[None] > 2.0:
+        # Both in high viscosity and high surface tension, for numerical stability it is recommend to set 0.0004
+        solver.dt[None] = ti.min(solver.dt[None], 0.0004)
+    gui.text('----------------------------')
+    gui.text('# of Fluid Particles')
+    gui.text('{}'.format(ps.total_fluid_particle_num))
+    gui.text('# of Rigid Particles')
+    gui.text('{}'.format(ps.total_rigid_particle_num))
+    gui.text('Total # of Particles')
+    gui.text('{}'.format(ps.total_particle_num))
     gui.end()
 
     scene.set_camera(camera)