# 复现iisph遇到问题

taichi 0.8.9

## 代码链接

1 Like

Hi @emt, 你好，非常欢迎你来到太极论坛。

1 Like

``````import taichi as ti
import math

poly6_factor = 315.0/(64.0*math.pi)

@ti.func
def Cubic(r_norm,h):
res = ti.cast(0.0, ti.f32)
k = 40 / 7 / math.pi/h**2
q = r_norm / h
if q <= 1.0:
if q <= 0.5:res = k * (6.0 * q**3 - 6.0 *q**2 + 1)
else:       res = k * 2 * ti.pow(1 - q, 3.0)
return res

@ti.func
res = ti.Vector([0.0,0.0])
k = 6. * 40 / 7 / math.pi / h ** 2
r_norm = r.norm()
q = r_norm / h
if r_norm > 1e-5 and q <= 1.0:
grad_q = r / (r_norm * h)
if q <= 0.5:     res = k * q * (3.0 * q - 2.0) * grad_q
else:            res = k * -(1.0 - q) ** 2 * grad_q
return res

dt  = 1e-4
r   = 0.05
h   = 4*r
V   = math.pi*r**2
rho0 = 1000
se_exp,se_coe = 7.0,50.0
fluid_rect = (20,120,1,1) #(width,hight,lowerleftX,lowerleftY)
bound      = (5,10)            # lowerleft bound is assumed (0,0)
n          = fluid_rect[0]*fluid_rect[1]
mass       = rho0*V
world2screen = 1.0/max(bound)

ti.init(arch = ti.cpu)

# particle properties
x   = ti.Vector.field(2,dtype = ti.f32 , shape = n)
xs  = ti.Vector.field(2,dtype = ti.f32 , shape = n) # screen position
v   = ti.Vector.field(2,dtype = ti.f32 , shape = n)
a   = ti.Vector.field(2,dtype = ti.f32 , shape = n)
rho = ti.field(dtype = ti.f32 , shape = n)
p   = ti.field(dtype = ti.f32 , shape = n)

# grid properties
cellNum,maxParticleInCell,maxNeighbor = 10007,100,100
gridData        = ti.field(dtype = ti.i32,shape = (cellNum,maxParticleInCell))
gridParticleNum = ti.field(dtype = ti.i32,shape = cellNum)
neighbor        = ti.field(dtype = ti.i32,shape = (n,maxNeighbor))
neighborNum     = ti.field(dtype = ti.i32,shape = n)

@ti.kernel
def Initialize():
for j in range(fluid_rect[1]):
for i in range(fluid_rect[0]):
v[j*fluid_rect[0]+i] = ti.Vector([0.0,0.0])
x[j*fluid_rect[0]+i] = ti.Vector([i*r,j*r])+ti.Vector([fluid_rect[2],fluid_rect[3]])

@ti.kernel
def StepBeforePressureSolve():
for i in gridParticleNum:gridParticleNum[i]=0
for i in x:
hashedIndex = (int(x[i][0]/2/h)*73856093+int(x[i][1]/2/h)*19349663)%cellNum
if gridParticleNum[hashedIndex]>=maxParticleInCell:continue
gridData[hashedIndex,idx] = i
for i in x:
neighborNum_i = 0
coord = int(x[i]/2/h)
for offset in ti.static([(-1,-1),(-1,0),(-1,1),(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)]):
nei_idx = ((coord[0]+offset[0])*73856093+(coord[1]+offset[1])*19349663)%cellNum
for celli in range(gridParticleNum[nei_idx]):
j = gridData[nei_idx,celli]
if i==j or (x[i]-x[j]).norm()>=h or neighborNum_i>=maxNeighbor:continue
neighbor[i,neighborNum_i]=j
neighborNum_i+=1
neighborNum[i] = neighborNum_i
for i in x:
rho_i = 0.0
for j in range(neighborNum[i]):
rho_i += mass*Cubic((x[i]-x[neighbor[i,j]]).norm(),h)
rho[i] = ti.max(rho_i,rho0)
for i in x:
ai    = ti.Vector([0.0,-9.81])
for nei_j in range(neighborNum[i]):
j = neighbor[i,nei_j]
a[i] = ai

#---------------------WCSPH-------------------------------------------------------------------------
@ti.kernel
def PressureSolve_WCSPH():
for i in x: p[i]  =  ti.max( se_coe*(ti.pow(rho[i]/rho0,se_exp)-1) ,0.0)

#----------------------IISPH-------------------------------------------------------------------------------
# A @ p = s
s     = ti.field(dtype = ti.f32,shape = n)
Adiag = ti.field(dtype = ti.f32,shape = n)
Ddiag  = ti.Vector.field(2,dtype = ti.f32,shape = n)
dijxpj = ti.Vector.field(2,dtype = ti.f32,shape = n)
p_new  = ti.field(dtype = ti.f32,shape = n)
rho_adv= ti.field(dtype = ti.f32,shape = n)

@ti.kernel
for i in v:
v[i] +=dt*a[i]
a[i]  = ti.Vector([0.0,0.0])
dii = ti.Vector([0.0,0.0])
for nei_j in range(neighborNum[i]):dii +=  CubicGradient(x[i]-x[neighbor[i,nei_j]],h)
Ddiag[i] = -dt*dt*mass/rho[i]**2*dii
for i in x:
aii = 0.0
# s_i = 0.0
for nei_j in range(neighborNum[i]):
j    = neighbor[i,nei_j]
aii+=(Ddiag[i] - dji).dot(dwij)
p[i] *= 0.5

@ti.kernel
def UpdatePressure_IISPH()->ti.f32: #return average error
for i in x:
dijpj = ti.Vector([0.0,0.0])
for nei_j in range(neighborNum[i]):
j    = neighbor[i,nei_j]
dijxpj[i] = dijpj
err = 0.0
# rho_avg = 0.0
for i in x:
sigma = 0.0
for nei_j in range(neighborNum[i]):
j = neighbor[i,nei_j]
sigma += mass*( dijxpj[i]-Ddiag[j]*p[j]-(dijxpj[j]- dji*p[i]) ).dot(dwij)
for i in p:p[i] = p_new[i]
ret = err/rho0/n
print(ret)
return ret

def PressureSolve_IISPH():
l = 0
while UpdatePressure_IISPH()>1e-3 and l<50: l+=1
#-----------------------------------------------------------------------------------------------------

@ti.kernel
def StepAfterPressureSolve():
for i in x: #compute pressure accelration
for nei_j in range(neighborNum[i]):
j = neighbor[i,nei_j]
for i in x:
v[i] += dt*a[i]
x[i] += dt*v[i]
for i in x:
if x[i][0]<r or x[i][0]>bound[0]-r:
x[i][0] = ti.max(ti.min(x[i][0],bound[0]-r),r)
v[i][0] = -0.9*v[i][0]
if x[i][1]<r or x[i][1]>bound[1]-r:
x[i][1] = ti.max(ti.min(x[i][1],bound[1]-r),r)
v[i][1] = -0.9*v[i][1]

@ti.kernel
for i in xs:xs[i] = x[i]*world2screen

def Main():
Initialize()
window = ti.ui.Window("SPH FLuid",(1000,1000))
canvas = window.get_canvas()
while not window.is_pressed(ti.ui.ESCAPE):
# for _ in range(1):
for _ in range(1):
StepBeforePressureSolve()
# PressureSolve_WCSPH()
PressureSolve_IISPH()
StepAfterPressureSolve()