# 关于SPH的理解

        density = [0 for _ in range(len(ps.pos))]
for idx, p in enumerate(ps.pos):
for nbIdx in ps.particle_neighbors[idx]:
pos = ps.pos[nbIdx]
dis = np.linalg.norm(np.array(p) - np.array(pos))


    def kernel(self, dis, h):

if self.dim == 2:
k = 40 / 7 / np.pi
else:
k = 8 / np.pi

q = dis/h
res = 0
if q <= 1.0:
if q <= 0.5:
q2 = q * q
q3 = q2 * q
res = k * (6.0 * q3 - 6.0 * q2 + 1)
else:
res = k * 2 *((1 - q)**3)
return res


import os.path

import numpy as np

class ParticleSystem:
# 用于记录粒子
def __init__(self):
self.dims = 2
self.pos = []  # 记录粒子位置
self.bounds = [] # minx, maxx, miny, maxy, minz, maxz
'''
记录粒子所在位置的压力，这个压力并不是被粒子带着走的
当粒子位置变化了之后，压力是需要更新的
'''
self.pressure = []
self.density = []
self.velocity = []

'''
在这个半径外，其他粒子对该粒子无影响
'''

'''
在计算密度，速度的微分等过程中，都需要计算邻居
'''
self.particle_neighbors = []

self.offset2 = [[x, y] for x in range(-1, 2) for y in range(-1, 2)]
self.offset3 = [[x, y, z] for x in range(-1, 2) for y in range(-1, 2) for z in range(-1, 2)]

def IsOutside(self, pos: list):
# 判断某个粒子是否在边界条件外
pass

def search_neighbors(self):
'''
但是，点太多了的话，求距离非常花时间
这里采用的方法，是先将区域画网格，然后在网格内的点才计算距离，减少计算
'''

gridNum = {}

for idx, p in enumerate(self.pos):
sub = self.pos_to_grid_sub(p)
index = self.sub_to_index(sub)
if index not in gridNum.keys():
gridNum[index] = []
gridNum[index].append(idx)

self.particle_neighbors = [[] for i in range(len(self.pos))]
for idx, p in enumerate(self.pos):
sub = self.pos_to_grid_sub(p)

if self.dims == 2:
offset = self.offset2
elif self.dims == 3:
offset = self.offset3

for of in offset:
newsub = [sub[0]+of[0], sub[1]+of[1]]
index = self.sub_to_index(newsub)
if index in gridNum.keys():
for ne in gridNum[index]:
distance = np.linalg.norm( np.array(p) - np.array(self.pos[ne]) )
if distance < self.support_radius and ne != idx:
self.particle_neighbors[idx].append(ne)

def pos_to_grid_sub(self, pos):

if not(( len(pos)==3 and len(self.bounds)==6 ) or (len(pos)==2 and len(self.bounds)==4)):
raise RuntimeError('coordinate wrong')

if self.dims == 2:
xmin, xmax, ymin, ymax = self.bounds

elif self.dims == 3:
xmin, xmax, ymin, ymax, zmin, zmax = self.bounds

gridX = int((pos[0] - xmin) / self.support_radius)
gridY = int((pos[1] - ymin) / self.support_radius)

if self.dims == 3:

if self.dims == 2:
return gridX, gridY
elif self.dims == 3:
return gridX, gridY, gridZ

def sub_to_index(self, sub):

if self.dims == 3:
index = sub[2]*self.Xlen*self.Ylen + sub[1]*self.Xlen + sub[0]
elif self.dims == 2:
index = sub[1] * self.Xlen + sub[0]

return index

class WaterParticleSystem(ParticleSystem):

def __init__(self):
super().__init__()
self.viscosity = 0.05  # viscosity
self.density0 = 1000

for x in np.arange(6, 9, 0.05):
for y in np.arange(2, 7, 0.05):
self.pos.append([x, y])
self.pressure.append(0) # 初始压力为0
self.density.append(self.density0)
self.velocity.append(0) # 初始速度
self.bounds = [0, 10, 0, 10]

def IsOutside(self, pos: list):
'''
定义二维边界，[0-10] x [0-10]
'''
x, y = pos
if x >= 0 and x <= 10 and y >=0 and y<=10:
return True
return False

class SPHSolver:

def __init__(self):
self.g = -0.98 # gravity
self.dt = 2e-4
self.dim = 2

def step(self, ps):
ps1 = self.substep(ps)
ps2 = self.enforce_boundary(ps1)
return ps2

def substep(self, ps):
self.compute_density(ps)

def enforce_boundary(self, ps):
pass

def compute_density(self, ps):

density = [0 for _ in range(len(ps.pos))]
for idx, p in enumerate(ps.pos):
for nbIdx in ps.particle_neighbors[idx]:
pos = ps.pos[nbIdx]
dis = np.linalg.norm(np.array(p) - np.array(pos))
print(np.max(density))

def compute_non_pressure_forces(self):
pass

def compute_pressure_forces(self):
pass

def kernel(self, dis, h):

if self.dim == 2:
k = 40 / 7 / np.pi
else:
k = 8 / np.pi

q = dis/h
res = 0
if q <= 1.0:
if q <= 0.5:
q2 = q * q
q3 = q2 * q
res = k * (6.0 * q3 - 6.0 * q2 + 1)
else:
res = k * 2 *((1 - q)**3)
return res

if __name__ == '__main__':

ps = WaterParticleSystem()
solver = SPHSolver()
ps.search_neighbors()
solver.compute_density(ps)



Python我也不了解，我试着看了一下，你的密度求解是用密度求和不是用连续方程，这部分没问题，但好像没有看到运动方程计算部分。sph方法是通过插值离散ns方程，ns方程是个偏微分方程，离散后可以求数值解

@zhang-qiang-github 这个cubic kernel的k是不是少除了一个 support_radius**self.dim?

    def kernel(self, dis, h):

if self.dim == 2:
k = 40 / 7 / np.pi
else:
k = 8 / np.pi

q = dis/h
res = 0
if q <= 1.0:
if q <= 0.5:
q2 = q * q
q3 = q2 * q
res = k * (6.0 * q3 - 6.0 * q2 + 1)
else:
res = k * 2 *((1 - q)**3)
return res


@mzhang 非常感谢细心的回答~~~

            for nbIdx in ps.particle_neighbors[idx]:
pos = ps.pos[nbIdx]
dis = np.linalg.norm(np.array(p) - np.array(pos))


import os.path

import numpy as np

class ParticleSystem:
# 用于记录粒子
def __init__(self):
self.dims = 2
self.pos = []  # 记录粒子位置
self.bounds = [] # minx, maxx, miny, maxy, minz, maxz
'''
记录粒子所在位置的压力，这个压力并不是被粒子带着走的
当粒子位置变化了之后，压力是需要更新的
'''
self.pressure = []
self.density = []
self.velocity = []

'''
在这个半径外，其他粒子对该粒子无影响
'''

'''
在计算密度，速度的微分等过程中，都需要计算邻居
'''
self.particle_neighbors = []

self.offset2 = [[x, y] for x in range(-1, 2) for y in range(-1, 2)]
self.offset3 = [[x, y, z] for x in range(-1, 2) for y in range(-1, 2) for z in range(-1, 2)]

def IsOutside(self, pos: list):
# 判断某个粒子是否在边界条件外
pass

def search_neighbors(self):
'''
但是，点太多了的话，求距离非常花时间
这里采用的方法，是先将区域画网格，然后在网格内的点才计算距离，减少计算
'''

gridNum = {}

for idx, p in enumerate(self.pos):
sub = self.pos_to_grid_sub(p)
index = self.sub_to_index(sub)
if index not in gridNum.keys():
gridNum[index] = []
gridNum[index].append(idx)

self.particle_neighbors = [[] for i in range(len(self.pos))]
for idx, p in enumerate(self.pos):
sub = self.pos_to_grid_sub(p)

if self.dims == 2:
offset = self.offset2
elif self.dims == 3:
offset = self.offset3

for of in offset:
newsub = [sub[0]+of[0], sub[1]+of[1]]
index = self.sub_to_index(newsub)
if index in gridNum.keys():
for ne in gridNum[index]:
distance = np.linalg.norm( np.array(p) - np.array(self.pos[ne]) )
if distance < self.support_radius and ne != idx:
self.particle_neighbors[idx].append(ne)

def pos_to_grid_sub(self, pos):

if not(( len(pos)==3 and len(self.bounds)==6 ) or (len(pos)==2 and len(self.bounds)==4)):
raise RuntimeError('coordinate wrong')

if self.dims == 2:
xmin, xmax, ymin, ymax = self.bounds

elif self.dims == 3:
xmin, xmax, ymin, ymax, zmin, zmax = self.bounds

gridX = int((pos[0] - xmin) / self.support_radius)
gridY = int((pos[1] - ymin) / self.support_radius)

if self.dims == 3:

if self.dims == 2:
return gridX, gridY
elif self.dims == 3:
return gridX, gridY, gridZ

def sub_to_index(self, sub):

if self.dims == 3:
index = sub[2]*self.Xlen*self.Ylen + sub[1]*self.Xlen + sub[0]
elif self.dims == 2:
index = sub[1] * self.Xlen + sub[0]

return index

class WaterParticleSystem(ParticleSystem):

def __init__(self):
super().__init__()
self.viscosity = 0.05  # viscosity
self.density0 = 1000

for x in np.arange(6, 9, 0.05):
for y in np.arange(2, 7, 0.05):
self.pos.append([x, y])
self.pressure.append(0) # 初始压力为0
self.density.append(self.density0)
self.velocity.append(0) # 初始速度
self.bounds = [0, 10, 0, 10]

def IsOutside(self, pos: list):
'''
定义二维边界，[0-10] x [0-10]
'''
x, y = pos
if x >= 0 and x <= 10 and y >=0 and y<=10:
return True
return False

class SPHSolver:

def __init__(self):
self.g = -0.98 # gravity
self.dt = 2e-4
self.dim = 2

def step(self, ps):
ps1 = self.substep(ps)
ps2 = self.enforce_boundary(ps1)
return ps2

def substep(self, ps):
self.compute_density(ps)

def enforce_boundary(self, ps):
pass

def compute_density(self, ps):

density = [0 for _ in range(len(ps.pos))]
for idx, p in enumerate(ps.pos):
for nbIdx in ps.particle_neighbors[idx]:
pos = ps.pos[nbIdx]
dis = np.linalg.norm(np.array(p) - np.array(pos))
print(np.max(density))

def compute_non_pressure_forces(self):
pass

def compute_pressure_forces(self):
pass

def kernel(self, dis, h):

if self.dim == 2:
k = 40 / 7 / np.pi
else:
k = 8 / np.pi
k /= h**self.dim

q = dis/h
res = 0
if q <= 1.0:
if q <= 0.5:
q2 = q * q
q3 = q2 * q
res = k * (6.0 * q3 - 6.0 * q2 + 1)
else:
res = k * 2 *((1 - q)**3)
return res

if __name__ == '__main__':

ps = WaterParticleSystem()
solver = SPHSolver()
ps.search_neighbors()
solver.compute_density(ps)


“密度求和”这种过程，我是不是可以理解成插值？

@mzhang 非常感谢提供的链接。

@mzhang 希望大侠能抽空解答一下困惑~~~

@RinSara @mzhang 给出的示例代码中，有一句代码是：self.m_V0 = 0.8 * self.particle_diameter ** self.dim，为什么计算质量是这么计算？这个0.8是个啥？

1. 为什么不用 \sum_j V_j W(r_i-r_j, h) 做归一化：

answer: 假如可以归一化，假设C有A/B两个邻居，那么公式为： f(C) = (\frac {V_A W(C-A, h)}{V_A W(C-A, h) + V_B W(C-B, h}f(A) + \frac {V_B W(C-B, h)}{V_A W(C-A, h) + V_B W(C-B, h}f(B)) 。也就是说，A对C的权重为 \frac {V_A W(C-A, h)}{V_A W(C-A, h) + V_B W(C-B, h)} ，C对A除的权重和却不是 V_A W(C-A, h) + V_B W(C-B, h) ，这会导致C对A的权重，与A对C的权重不一样。这不符合物理原理

1. 为什么可以用 r_j 来插值 r_i:

answer:t_n 时刻，已有位置 p_{i=1, 2, ..., m}^{t_n} 的密度。然后在 t_{n+1} ，需要计算 m 个粒子的密度。由于在 t_{n+1} 时刻，粒子的位置已经移动到了: p_{i=1, 2, ..., m}^{t_{n+1}} 。为了计算 p_i^{t_{n+1}} ，用其周围的 p_j^{t_n} 来计算。

1 Like

（二维中） 0.8
（2r）^2=（pi*r^2）