大家好,最近在使用Taichi上面測試demo的mpm99物質點法程式碼出了點問題想請教大神們
我試圖在原始的mpm99上面去修改kernel function以及support size,但在修改途中發現下面兩個方式得到的結果完全不一樣:
在mpm99中grid to particle (G2P) 過程中我根據網上查到的kernel的原始定義定義添加了額外的for-loop來計算
w_old[0][0] = 0.5 * (1.5 - fx[0])**2
w_old[0][1] = 0.5 * (1.5 - fx[1])**2
w_old[1][0] = 0.75 - (fx[0] - 1)**2
w_old[1][1] = 0.75 - (fx[1] - 1)**2
w_old[2][0] = 0.5 * (fx[0] - 0.5)**2
w_old[2][1] = 0.5 * (fx[1] - 0.5)**2
for idx, idy in ti.static(ti.ndrange(numNeighbor, dim)): # Loop over possible neighbor in x y direction
distReg = abs(fx[idy]-float(idx))# this is the x or y in the kernel expression, the reason for 1 is because of the float of idx
if (distReg) >= 0 and (distReg) <= 0.5:
w[idx][idy] = 0.75-distReg**2
elif (distReg) > 0.5 and (distReg) <= 1.5:
w[idx][idy] = 0.5*(1.5-distReg)**2
elif (distReg) > 1.5:
w[idx][idy] = 0
其中 w_old 是原本而mpm99提供的Quadratic kernels,我試圖透過上面的方式去更改使用不同的kernel,但跑出來的結果卻是非常不穩定很快便發散掉了(程式執行不會報錯)。
因此我試圖在for-loop裡面增加一個if判斷來檢查是哪一個kernel的數學式子可能寫錯,卻發現結果能夠順利運行且數值穩定,得到跟原始的mpm99一致的結果(code如下)
w_old[0][0] = 0.5 * (1.5 - fx[0])**2
w_old[0][1] = 0.5 * (1.5 - fx[1])**2
w_old[1][0] = 0.75 - (fx[0] - 1)**2
w_old[1][1] = 0.75 - (fx[1] - 1)**2
w_old[2][0] = 0.5 * (fx[0] - 0.5)**2
w_old[2][1] = 0.5 * (fx[1] - 0.5)**2
for idx, idy in ti.static(ti.ndrange(numNeighbor, dim)): # Loop over possible neighbor in x y direction
distReg = abs(fx[idy]-float(idx))# this is the x or y in the kernel expression, the reason for 1 is because of the float of idx
if (distReg) >= 0 and (distReg) <= 0.5:
w[idx][idy] = 0.75-distReg**2
elif (distReg) > 0.5 and (distReg) <= 1.5:
w[idx][idy] = 0.5*(1.5-distReg)**2
elif (distReg) > 1.5:
w[idx][idy] = 0
# double check if the proposed algorithm matches the original setting, if not, output the information**
if abs(w_old[idx][idy]-w[idx][idy]) > 1e-4: **
print('x[p] :',x[p])**
print('base :',base)**
print('distReg :',distReg)**
print('fx[0] :',fx[idy])**
print('w_old[idx][idy]',w_old[idx][idy])**
print('w[idx][idy]',w[idx][idy])**
想請教大神們在taichi裡面巢狀for-loop 再額外添加for-loop或if結構有特別的語法限制嗎?我是用單個thread下去執行,但不管單個或多個thread得到的結果都是一樣的結論。我看了網上的視頻只記得會影響併行計算。