在实现混合法时,p2g,g2p用了pic ,用CG求解了压力泊松方程。
运行时会爆炸,debug时出现错误:
出错的代码块;
@ti.kernel
def ParticlesToGrid():
#分别转化 U,V。先尝试用Quadratic B-spline进行插值
for p in x:
base_u = (x[p] * inv_dx -ti.Vector([0.0,0.5]) ).cast(int)
base_v=(x[p] * inv_dx -ti.Vector([0.5,0.0]) ).cast(int)
fx_u = x[p] * inv_dx - base_u.cast(float)
fx_v = x[p] * inv_dx - base_v.cast(float)
# Quadratic B-spline
w_u = [0.5 * (1.5 - fx_u) ** 2, 0.75 - (fx_u - 1) ** 2, 0.5 * (fx_u - 0.5) ** 2]
w_v = [0.5 * (1.5 - fx_v) ** 2, 0.75 - (fx_v - 1) ** 2, 0.5 * (fx_v - 0.5) ** 2]
# 分别处理 U,V
for i in ti.static(range(3)):
for j in ti.static(range(3)):
weight_u = w_u[i][0] * w_u[j][1]
weight_v = w_v[i][0] * w_v[j][1]
index_u = base_u + ti.Vector([i, j])
index_v = base_v + ti.Vector([i, j])
if index_u[0] <= n_grid and index_u[1] <= n_grid - 1:
grid_u[index_u] += weight_u * v[p][0]
grid_m_u[index_u] += weight_u
if index_v[0] <= n_grid - 1 and index_v[1] <= n_grid:
grid_v[index_v] += weight_v * v[p][1]
grid_m_v[index_v] += weight_v
for i, j in grid_m_u:
if grid_m_u[i, j] > 0:
inv_m = 1 / grid_m_u[i, j]
grid_u[i, j] = inv_m * grid_u[i, j]
for i, j in grid_m_v:
if grid_m_v[i, j] > 0:
inv_v = 1 / grid_m_v[i, j]
grid_v[i, j] = inv_v * grid_v[i, j]
就是这个u方向的速度,P2G时出现了问题,怎么也改不对,求大佬指导