关于APIC仿真水平线条掉落中精度的疑问

想请教各位,我用 MPM + APIC 方法仿真简单的顶点掉落,发现掉落高度不一样。
场景只有一根水平放置的线条,顶点均匀采样生成。我发现掉落过程中顶点掉落的高度会有一些微小的偏差,这是浮点数误差的问题吗?还是有可能在原理层面解决掉的?

附上代码,我计算相邻顶点的高度差,每 100 帧 log 一次。开始 200 帧没出现高度不一致,继续仿真出现偏差。

# most code from https://zhuanlan.zhihu.com/p/414356129

import taichi as ti
import random
import math

ti.init(arch = ti.cpu)

dim = 2
n_grid = 256
dx = 1 / n_grid
inv_dx = 1 / dx
dt = 4.0e-5

n_type2 = 200
Length = 0.75
sl = Length/ (n_type2-1)
volume2 =  dx*Length / (n_type2)

p_rho = 1

start_pos = ti.Vector([0.2, 0.8])

x2 = ti.Vector.field(2, dtype=float, shape=n_type2) # position
v2 = ti.Vector.field(2, dtype=float, shape=n_type2) # velocity
C2 = ti.Matrix.field(2, 2, dtype=float, shape=n_type2) # affine velocity field

grid_v = ti.Vector.field(2, dtype= float, shape=(n_grid, n_grid))
grid_m = ti.field(dtype=float, shape=(n_grid, n_grid))

@ti.kernel
def Reset():
    for i, j in grid_m:
        grid_v[i, j] = [0, 0]
        grid_m[i, j] = 0

@ti.kernel
def Particle_To_Grid():
    for p in x2:
        base = (x2[p] * inv_dx - 0.5).cast(int)
        fx = x2[p] * inv_dx - base.cast(float)
        w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5) ** 2]
        affine = C2[p]
        mass = volume2* p_rho
        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i,j])
            weight = w[i][0]*w[j][1]
            grid_m[base + offset] += weight * mass
            dpos = (offset.cast(float) - fx) * dx
            grid_v[base + offset] += weight * mass * (v2[p] +  affine@dpos)

bound = 3
@ti.kernel
def Grid_Collision():
    for i, j in grid_m:
        if grid_m[i, j] > 0:
            grid_v[i, j] /= grid_m[i, j]
            grid_v[i, j].y -= dt * 9.80

            if i < bound and grid_v[i, j].x < 0:
                grid_v[i, j].x = 0
            if i > n_grid - bound and grid_v[i, j].x > 0:
                grid_v[i, j].x = 0
            if j < bound and grid_v[i, j].y < 0:
                grid_v[i, j].y = 0
            if j > n_grid - bound and grid_v[i, j].y > 0:
                grid_v[i, j].y = 0

@ti.kernel
def Grid_To_Particle():
    for p in x2:
        base = (x2[p] * inv_dx - 0.5).cast(int)
        fx = x2[p] * inv_dx - base.cast(float)
        w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
        new_v = ti.Vector.zero(float, 2)
        new_C = ti.Matrix.zero(float, 2, 2)
        new_x = ti.Vector.zero(float, 2)
        for i, j in ti.static(ti.ndrange(3, 3)):
            dpos = ti.Vector([i, j]).cast(float) - fx
            g_v = grid_v[base + ti.Vector([i, j])]
            weight = w[i][0] * w[j][1]
            new_v += weight * g_v
            new_C += 4 * weight * g_v.outer_product(dpos) * inv_dx
            offset = ti.Vector([i, j]).cast(float)
            new_x += weight * dx * (base.cast(float) + offset)
        v2[p] = new_v
        x2[p] += dt * v2[p]
        C2[p] = new_C

@ti.kernel
def initialize():
    for i in range(n_type2):
        x2[i] = ti.Vector([start_pos[0]+ i * sl, start_pos[1]])
        v2[i] = ti.Matrix([0, 0])
        C2[i] =  ti.Matrix([[0,0],[0,0]])

def main():
    initialize()
    frame = 0
    for __ in range(100):
        maximal = 0.
        for _ in range(100):
            Reset()
            Particle_To_Grid()
            Grid_Collision()
            Grid_To_Particle()
            x2_np = x2.to_numpy()
            for i in range(n_type2 - 1):
                d = x2_np[i+1] - x2_np[i]
                maximal = max(maximal, math.fabs(d[1]))
            frame += 1

        print("max difference " + str(maximal) + "\tduring frame[" + str(frame - 100)+ "," +str(frame) + "]", flush=True)

if __name__ == "__main__":
    main()

运行结果

[Taichi] version 0.8.4, llvm 10.0.0, commit 895881b5, linux, python 3.9.9
[Taichi] Starting on arch=x64
max difference 0.0	during frame[0,100]
max difference 0.0	during frame[100,200]
max difference 5.960464477539063e-08	during frame[200,300]
max difference 5.960464477539063e-08	during frame[300,400]
max difference 5.960464477539063e-08	during frame[400,500]
max difference 5.960464477539063e-08	during frame[500,600]
max difference 5.960464477539063e-08	during frame[600,700]
max difference 5.960464477539063e-08	during frame[700,800]
max difference 1.1920928955078125e-07	during frame[800,900]
max difference 1.1920928955078125e-07	during frame[900,1000]
max difference 1.1920928955078125e-07	during frame[1000,1100]
max difference 1.1920928955078125e-07	during frame[1100,1200]
max difference 1.1920928955078125e-07	during frame[1200,1300]
max difference 1.1920928955078125e-07	during frame[1300,1400]
max difference 2.384185791015625e-07	during frame[1400,1500]
max difference 2.384185791015625e-07	during frame[1500,1600]
max difference 2.384185791015625e-07	during frame[1600,1700]
max difference 2.384185791015625e-07	during frame[1700,1800]
max difference 2.384185791015625e-07	during frame[1800,1900]
max difference 2.384185791015625e-07	during frame[1900,2000]
max difference 2.980232238769531e-07	during frame[2000,2100]
max difference 4.172325134277344e-07	during frame[2100,2200]
max difference 4.172325134277344e-07	during frame[2200,2300]
max difference 4.172325134277344e-07	during frame[2300,2400]
max difference 4.172325134277344e-07	during frame[2400,2500]
max difference 4.172325134277344e-07	during frame[2500,2600]
max difference 5.960464477539062e-07	during frame[2600,2700]
max difference 6.556510925292969e-07	during frame[2700,2800]
max difference 6.556510925292969e-07	during frame[2800,2900]
max difference 6.556510925292969e-07	during frame[2900,3000]
max difference 7.152557373046875e-07	during frame[3000,3100]
max difference 9.5367431640625e-07	during frame[3100,3200]
max difference 9.5367431640625e-07	during frame[3200,3300]
max difference 9.5367431640625e-07	during frame[3300,3400]
max difference 9.5367431640625e-07	during frame[3400,3500]
max difference 1.1920928955078125e-06	during frame[3500,3600]
max difference 1.3113021850585938e-06	during frame[3600,3700]
max difference 1.3113021850585938e-06	during frame[3700,3800]
max difference 1.3113021850585938e-06	during frame[3800,3900]
max difference 1.430511474609375e-06	during frame[3900,4000]
max difference 1.7285346984863281e-06	during frame[4000,4100]
max difference 1.7285346984863281e-06	during frame[4100,4200]
max difference 1.7285346984863281e-06	during frame[4200,4300]
max difference 1.7881393432617188e-06	during frame[4300,4400]
max difference 2.1457672119140625e-06	during frame[4400,4500]
max difference 2.2649765014648438e-06	during frame[4500,4600]
max difference 2.2649765014648438e-06	during frame[4600,4700]
max difference 2.2649765014648438e-06	during frame[4700,4800]
max difference 2.5033950805664062e-06	during frame[4800,4900]
max difference 2.86102294921875e-06	during frame[4900,5000]
max difference 2.86102294921875e-06	during frame[5000,5100]
max difference 2.86102294921875e-06	during frame[5100,5200]
max difference 3.039836883544922e-06	during frame[5200,5300]
max difference 3.3974647521972656e-06	during frame[5300,5400]
max difference 3.4570693969726562e-06	during frame[5400,5500]
max difference 3.4570693969726562e-06	during frame[5500,5600]
max difference 3.5762786865234375e-06	during frame[5600,5700]
max difference 3.874301910400391e-06	during frame[5700,5800]
max difference 4.231929779052734e-06	during frame[5800,5900]
max difference 4.231929779052734e-06	during frame[5900,6000]
max difference 4.291534423828125e-06	during frame[6000,6100]
max difference 4.589557647705078e-06	during frame[6100,6200]
max difference 4.827976226806641e-06	during frame[6200,6300]
max difference 5.0067901611328125e-06	during frame[6300,6400]
max difference 5.27501106262207e-06	during frame[6400,6500]
max difference 5.424022674560547e-06	during frame[6500,6600]
max difference 5.692243576049805e-06	during frame[6600,6700]
max difference 5.900859832763672e-06	during frame[6700,6800]
max difference 6.198883056640625e-06	during frame[6800,6900]
max difference 6.407499313354492e-06	during frame[6900,7000]
max difference 6.645917892456055e-06	during frame[7000,7100]
max difference 6.9141387939453125e-06	during frame[7100,7200]
max difference 7.152557373046875e-06	during frame[7200,7300]
max difference 7.450580596923828e-06	during frame[7300,7400]
max difference 7.659196853637695e-06	during frame[7400,7500]
max difference 7.957220077514648e-06	during frame[7500,7600]
max difference 8.165836334228516e-06	during frame[7600,7700]
max difference 8.52346420288086e-06	during frame[7700,7800]
max difference 8.791685104370117e-06	during frame[7800,7900]
max difference 9.059906005859375e-06	during frame[7900,8000]
max difference 9.357929229736328e-06	during frame[8000,8100]
max difference 9.685754776000977e-06	during frame[8100,8200]
max difference 1.0013580322265625e-05	during frame[8200,8300]
max difference 1.0371208190917969e-05	during frame[8300,8400]
max difference 1.0713934898376465e-05	during frame[8400,8500]
max difference 1.1026859283447266e-05	during frame[8500,8600]
max difference 1.1354684829711914e-05	during frame[8600,8700]
max difference 1.1667609214782715e-05	during frame[8700,8800]
max difference 1.2040138244628906e-05	during frame[8800,8900]
max difference 1.2382864952087402e-05	during frame[8900,9000]
max difference 1.271069049835205e-05	during frame[9000,9100]
max difference 1.30385160446167e-05	during frame[9100,9200]
max difference 1.3396143913269043e-05	during frame[9200,9300]
max difference 1.3738870620727539e-05	during frame[9300,9400]
max difference 1.411139965057373e-05	during frame[9400,9500]
max difference 1.446157693862915e-05	during frame[9500,9600]
max difference 1.4819204807281494e-05	during frame[9600,9700]
max difference 1.519918441772461e-05	during frame[9700,9800]
max difference 1.5575438737869263e-05	during frame[9800,9900]
max difference 1.5966594219207764e-05	during frame[9900,10000]