Mpm99, for-loop 裡面添加 if 判斷式造成結果不一致

大家好,最近在使用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得到的結果都是一樣的結論。我看了網上的視頻只記得會影響併行計算。

你好,欢迎来到Taichi社区!nested for loop的影响一般来说是inner for loop会被串行执行,仿真结果不同的话只看这部分代码有点难以定位问题,这边方便附上完整代码供调试吗? :smiley: