想请教各位，我用 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]