运行环境
taichi 0.8.9
问题描述
我在课程提供的sph_base和particle_system基础上，复现iisph原论文时，发现无法收敛压强，具体表现为a_ii的值很小，迭代过程中压强p_i,l+1 = (density_0 - desity_i,adv - (…))/a_ii非常大，超过正常的压强应有的数值
我读公式没发现问题，求大佬们教教
1 个赞
不好意思，我之前疏忽忘记放了 ……
已经更新在repo里了，运行ii_demo.py会调用iisph.py和其他文件，麻烦大佬看一下谢谢
1 个赞
最后问题解决了吗？我最近写iisph时也出现了不收敛的问题。
import taichi as ti
import math
poly6_factor = 315.0/(64.0*math.pi)
spiky_grad_factor = -45.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
def CubicGradient(r,h):
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
idx = ti.atomic_add(gridParticleNum[hashedIndex],1)
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]
ai += V*(v[i]-v[j]).dot(x[i]-x[j])/((x[i]-x[j]).norm_sqr()+0.01*h*h)*CubicGradient(x[i]-x[j],h)
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
def FillADiagonalAnds_IISPH():
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
rho_adv[i] = rho[i]
for nei_j in range(neighborNum[i]):
j = neighbor[i,nei_j]
dwij = CubicGradient(x[i]-x[j],h)
dji = -dt**2*mass*CubicGradient(x[j]-x[i],h)/rho[i]**2
aii+=(Ddiag[i] - dji).dot(dwij)
rho_adv[i] += dt*mass*(v[i]-v[j]).dot(dwij)
Adiag[i] = aii*mass
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]
dijpj += (-dt*dt*mass/rho[j]**2)*CubicGradient(x[i]-x[j],h)*p[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]
dwij = CubicGradient(x[i]-x[j],h)
dji = -dt**2*mass*CubicGradient(x[j]-x[i],h)/rho[i]**2
sigma += mass*( dijxpj[i]-Ddiag[j]*p[j]-(dijxpj[j]- dji*p[i]) ).dot(dwij)
p_new[i] = ti.max( 0.5*p[i]+0.5*(rho0-rho_adv[i]-sigma)/Adiag[i],0.0)
err += (sigma + Adiag[i]*p_new[i] -(rho0-rho_adv[i]))
for i in p:p[i] = p_new[i]
ret = err/rho0/n
print(ret)
return ret
def PressureSolve_IISPH():
l = 0
FillADiagonalAnds_IISPH()
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]
a[i] += -mass*(p[j]/rho[j]**2 + p[i]/rho[i]**2)*CubicGradient(x[i]-x[j],h)
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
def UpdateScreenPosition():
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()
UpdateScreenPosition()
canvas.set_background_color((245/255,245/255,245/255))
canvas.circles(xs,radius = r*world2screen/2,color=(183/255,134/255,11/255))
window.show()
Main()