[bug] A possible hidden bug of MPM examples

[chunleili]

Describe the bug
In taichi MPM official examples, the quadratic interpolation weight should be non-negative. But a zero position particle will get a negative result. A possible reason is that base = ti.cast(Xp - 0.5, ti.i32) will get 0, causing fx = Xp - base =0 and making 0.75 - (fx - 1)**2<0.

This bug will not cause error in MPM simulation because there is a separate discussion of the edge case later. But if we run the p2g separately, it will produce the bug. This bug may cause troubles for those who copy the official examples and making it as reference. Therefore it is worthwhile to fix this bug.

To Reproduce
Please post a minimal sample code to reproduce the bug.
The developer team will put a higher priority on bugs that can be reproduced within 20 lines of code. If you want a prompt reply, please keep the sample code short and representative.

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()

Log/Screenshots

weight<0! [0, 0, 0]
weight<0! [0, 0, 0]
weight<0! [0, 0, 0]
...

Possible solution
change to
base = ti.cast(ti.floor(Xp - 0.5), ti.i32)

[lin-hitonami]

Thanks! @chunleili Would you like to open a PR to fix that?

[yuanming-hu]

You are right - we should use ti.floor(..., ti.i32).