已挂issue:
[bug] A possible hidden bug of MPM examples · Issue #7792 · taichi-dev/taichi (github.com)
描述错误
在太极MPM官方示例中，二次插值权重应为非负值。但是零位置粒子会得到负结果。一个可能的原因是base = ti.cast(Xp - 0.5, ti.i32)将得到 0，导致fx = Xp - base =0 , 从而0.75 - (fx - 1)**2<0
此错误不会在 MPM 模拟中导致错误，因为稍后会单独讨论边缘情况。但是如果我们单独运行 p2g，它会产生错误。此错误可能会给那些复制官方示例并将其作为参考的人带来麻烦。
BUG复现
import taichi as ti
ti.init(arch=ti.gpu)
n_particles = 100
dim, n_grid = 3, 32
x = ti.Vector.field(dim, float, n_particles)
grid_m = ti.field(float, (n_grid, ) * dim)
neighbour = (3, ) * dim
@ti.kernel
def test():
dx = 0.1
p_mass = 1.0
for p in x:
Xp = x[p] / dx
base = ti.cast(Xp - 0.5, ti.i32) #bug at 0,0,0
# base = ti.cast(ti.floor(Xp - 0.5), ti.i32) #possible fix of the bug
fx = Xp - base
w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
for offset in ti.static(ti.grouped(ti.ndrange(*neighbour))):
weight = 1.0
for i in ti.static(range(dim)):
weight *= w[offset[i]][i]
grid_m[base + offset] += weight * p_mass
if(weight < 0):
print('weight<0!', base)
if __name__ == '__main__':
test()
output
weight<0! [0, 0, 0]
weight<0! [0, 0, 0]
weight<0! [0, 0, 0]
可能的解决方案
改为
base = ti.cast(ti.floor(Xp - 0.5), ti.i32)
