请问关于mpm88的fx和base到底代表什么?

各位好
关于mpm88里面,总是先有这么一行
base=int(x[p]/dx-0.5)
fx=x[p]/dx-base
请教一下这两个变量到底是什么含义啊?

为什么不直接用x而要非得减去0.5呢?

而且后面的形函数w,为什么要带入X=fx, fx-1, fx-2啊?

后面又为什么非得用w[0]*w[1]?

3 个赞

您好,首先非常欢迎您加入Taichi社区。

关于您的问题,这涉及到MPM的具体数学公式。

  1. 关于base和fx的求解,可以参考here 的公式121.
  2. 其次形函数w的使用,可以参考上面文档的Page 33.
1 个赞

您好
感谢您的回答

我认真读过您所发的 蒋老师的 MPM course
然而
我所疑问的是
为什么该插值在代码中是这样处理的。关于这一点,我思考了很多天都没有理解。
尤其是base 和fx的具体几何含义。如果可以,可否用图片标注fx和base的具体所指?

以及C++注释版本中,提到把fx fx-1 fx-2分别带入公式的x中,我实在不能理解为什么这样做。

祝好

嗯嗯,这是一个很好的问题,确实很容易让你迷惑。
首先,减去0.5的原因是我们一般认为粒子在网格的中间。胡老师在Games 201, Lec 7, page 16,讲过为什么减去0.5,可以在B站具体听一下讲解here

MPM中P2G阶段,需要将particle上的信息按照权重 w ,分散到周围的网格点上,这个例子里面就是 3\times3 的网格大小。
base 就是 这个网格的左下角网格点的位置。fx就是particle 到base网格点的差值。

我们可以以 p=[1.2,0.6], dx = 1.0 为例,研究一下权重 w 的求解。

我们这个例子中使用的是Quadratic kernels

那么 px 轴上,根据距离差可以算出和对应网格顶点的权重
base的影响(权重 w )应该是 N(1.2)
B的影响(权重 w )应该是 N(0.2)
C的影响(权重 w )应该是 N(-0.8)

你现在应该都还是很清楚的,但是你会看到这个公式和代码中的公式是不一样的。这就要涉及到一些数学了。我们将核函数 N 根据定义域分为三个部分 [-1.5, -0.5], [-0.5, 0.5], [0.5,1.5] 。这三个部分都可以用 [0.5,1.5] 的三个函数来表示。依然用图片给你展示一下:

  1. N(x) \in [0.5,1.5] = f(x) \in [0.5,1.5]

  2. N(x) \in [-0.5,0.5] = g(x) \in [0.5,1.5]

  3. N(x) \in [-1.5,-0.5] = h(x) \in [0.5,1.5]

综上,你可以看出来
base的影响(权重 w )应该是 N(1.2), 也就是 f(1.2)
B的影响(权重 w )应该是 N(0.2), 也就是 g(1.2)
C的影响(权重 w )应该是 N(-0.8) , 也就是 h(1.2)
其中 1.2 = fx.x
那么 mpm88.py的Line 37

w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]

中的 w[0].x = f(1.2) = N(1.2), w[1].x = g(1.2) = N(0.2) , w[2].x = h(1.2) = N(-0.8)

Y轴同理 w[0].y = f(0.6)=N(0.6), w[1].y = g(0.6)=N(-0.4), w[2].y = h(0.6) = N(-1.4)

7 个赞

过这么久了,自己动手写了一下SPH,我总算是明白为什么会减去0.5了。

不知道我想的对不对,如果不对请指正。

这是因为粒子在右边界和上边界计算出的网格编号会溢出,减去0.5相当于把整张网格向左下方平移了0.5个单位。

为什么会溢出,是这样的:
假设整个计算区域大小是100x100, 网格尺寸是1x1

我们利用int()来向下取整,假如粒子在上边界,比如position=(50,100),那么,当粒子的y坐标100转换为网格ID的时候,由于是向下取整的,网格的y编号就会变为100,但是网格的y编号的范围是0~99(因为从0开始算100个网格)。这就导致计算的网格编号超出范围了!

有边界也是同理的。

假如不是恰好在边界,而是稍微超过边界,比如100.1,那也会造成溢出。

但是下边界和左边界就没事,因为-0.1取整之后还是0,没有超出范围。所以左下边界是本身就具有一个单位的裕量的。

解决方案就是把整个网格整体向左下平移0.5个单位。这样原本100.1也会变成99.6,取整之后变为99,不会超出范围。移动之后,上下左右四个边界都拥有了半个单位的裕量。

同时,还应该注意边界处理的时候,把超过边界的值挪动到离着边界还有一个裕量的位置,而不是恰好在边界上,比如100-eps,eps可以设置为0.1

同时,最好给边界多预留一些距离,比如MPM99里面处理边界的时候,就在离着边界的距离小于3个单位的时候,就假设它穿过了边界,然后处理。

3 个赞

感谢解答,对于理解代码很有用。

当0.5<=x<1.5, base node 为0
当1.5<=1.5<2.5, base node 为1

但是我有一个疑问, 当粒子 x<0.5时,base node 依然为0,但是weight的计算(w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2])并不适用,对于这些粒子,怎么处理?

Hi @xiangzi , 应该还是适用的。你可以跟着代码和公式过一遍。即使base node 还是0, 也还是可以以base node 为原点,遍历3x3的网格顶点,而对应的weight计算公式还是不变。

假设找到一个粒子,使得fx=0.1, base node 是0, 那么把fx代入w计算得到的weight是=[0.98, -0.06, 0.08], 这是不合理的,是我哪里理解错了吗

请问,计算出来的fx在x方向是fx.x=1.2,然后带入w;但为什么y方向的fx.y却是不一样的值f(0.6),g(1.2),h(1.2);以x的f(x),g(x),h(x)来计算,在计算出来的y方向的不应该是fx.y=0.6,f(0.6),g(0.6),h(0.6)这样吗?

1 个赞

系统梳理了一遍原理~如有不严谨的地方请多指教 :grin:

/(ㄒoㄒ)/~~不让俺上传多张图片,那就只能传到B站上了

https://www.bilibili.com/read/cv16719283

1 个赞

原理与程序部分:

例子:
https://www.bilibili.com/read/cv16719283

2 个赞

同学,谢谢你。我仔细看看!

谢谢指出错误,已经更改了。

大佬您好,我带了几个点检验了一下胡老师代码(pic_vs_apic.py)里计算base的方法,得到的结果似乎是B点(借用一下您的图)


当dx等于网格维度的倒数的时候,好像base是当前位置左下角的网格点。这里有点没搞明白。。。
于是我按照自己的理解写了一个计算base的方法。base是距离该位置最近的网格点,然后以base为中心计算9个点的权重,最后完成P2G和G2P操作。如图所示,红框是以D点为中心的网格,在红框范围内的粒子计算出来的base就是D点。

代码如下:

@ti.kernel
def my_PIC():
    for p in x:
        # 胡老师的算法
        # base = (x[p] * inv_dx - 0.5).cast(int)
        # fx = x[p] * inv_dx - base.cast(float)
        # Quadratic B-spline
        # w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
        # 我的理解
        base = ti.round((x[p] * inv_dx - 0.5), int)  # 根据粒子位置计算最近网格位置,以该位置作为中心网格取3*3个网格点分别计算权值
        fx = x[p] * inv_dx - base.cast(float) - 0.5  # 粒子位置-中心网格位置=粒子偏移中心网格的量
        # 以x轴为例,左中右与粒子位置的距离分别为:左=1+fx.x , 中=fx.x , 右=1-fx.x
        # 又因为-0.5<fx.x<0.5,所以0.5<1+fx.x<1.5,0.5<1-fx.x<1.5
        # 所以w计算如下
        w = [0.5 * (0.5 - fx) ** 2, 0.75 - fx ** 2, 0.5 * (0.5 + fx) ** 2]
        for i in ti.static(range(3)):
            for j in ti.static(range(3)):
                # offset = ti.Vector([i, j]) # 胡老师的算法
                offset = ti.Vector([i - 1, j - 1])  # 我的理解:因为base是中心点,所以他的左下邻居的偏移量应该是[-1,-1]
                weight = w[i][0] * w[j][1]
                grid_v[base + offset] += weight * v[p]
                grid_m[base + offset] += weight

    for i, j in grid_m:
        if grid_m[i, j] > 0:
            inv_m = 1 / grid_m[i, j]
            grid_v[i, j] = inv_m * grid_v[i, j]

    for p in x:
        # base = (x[p] * inv_dx - 0.5).cast(int)
        # fx = x[p] * inv_dx - base.cast(float)
        # Quadratic B-spline
        # w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
        # 我的理解
        base = ti.round((x[p] * inv_dx - 0.5), int)  # 根据粒子位置计算最近网格位置,以该位置作为中心网格取3*3个网格点分别计算权值
        fx = x[p] * inv_dx - base.cast(float) - 0.5  # 粒子位置-网格中心位置=粒子偏移中心网格的量
        w = [0.5 * (0.5 - fx) ** 2, 0.75 - fx ** 2, 0.5 * (0.5 + fx) ** 2]
        new_v = ti.Vector.zero(ti.f32, 2)
        for i in ti.static(range(3)):
            for j in ti.static(range(3)):
                weight = w[i][0] * w[j][1]
                # new_v += weight * grid_v[base + ti.Vector([i, j])] # 胡老师的算法
                new_v += weight * grid_v[base + ti.Vector([i - 1, j - 1])]  # 我的理解:因为base是中心点,所以他的左下邻居的偏移量应该是[-1,-1]

        x[p] = clamp_pos(x[p] + v[p] * dt)
        v[p] = new_v

效果看起来和胡老师的代码是一样的。
我有两个疑问:
1.我的理解到底对不对。
2.如果我的理解是对的,那么我的方法和胡老师的方法在效率上哪个更好一点。
感谢大佬解惑!

1 个赞

base 点的坐标,你可以在python里跑一下看看

base.x = int(1.2 - 0.5) = int(0.7)  = 0
base.y = int(0.6 - 0.5) = int(0.1)  = 0

减0.5和的问题我也困扰了好久,貌似想通了,做了下记录。
不知道我这个理解对不对。。。求证实







10 个赞

非常详细的解释哈,确实如你所说,是在黑色格点上计算权重的。

您好 请问 在flip 中为何一般都是在x y z 三个方向上为何是用 4X3X3 而不是 3X3X3来算权重?

@ti.func
def interp_grid(base, frac, vp):
    
    # Index on sides
    #idx_side = ti.Vector([base-1, base, base+1, base+2])
    idx_side = [base-1, base, base+1, base+2]

    # Weight on sides
    #w_side = ti.Vector([quadratic_kernel(1.0+frac), quadratic_kernel(frac), quadratic_kernel(1.0-frac), quadratic_kernel(2.0-frac)])
    w_side = [quadratic_kernel(1.0+frac), quadratic_kernel(frac), quadratic_kernel(1.0-frac), quadratic_kernel(2.0-frac)]
    
    # Index on center 
    #idx_center = ti.Vector([base-1, base, base+1])
    idx_center = [base-1, base, base+1]

    # weight on center 
    #w_center= ti.Vector([quadratic_kernel(0.5+frac), quadratic_kernel(ti.abs(0.5-frac)), quadratic_kernel(ti.abs(0.5-frac)), quadratic_kernel(1.5-frac)])
    w_center= [quadratic_kernel(0.5+frac), quadratic_kernel(ti.abs(0.5-frac)), quadratic_kernel(1.5-frac)]

    #for i, j, k in ti.ndrange(4, 3, 3):
    for i in ti.static(range(4)):
        for j in ti.static(range(3)):
            for k in ti.static(range(3)):
                w = w_side[i].x * w_center[j].y * w_center[k].z
                idx = (idx_side[i].x, idx_center[j].y, idx_center[k].z)
                grid_v_x[idx] += vp.x * w
                grid_w_x[idx] += w
        

    #for i, j, k in ti.ndrange(3, 4, 3):
    for i in ti.static(range(3)):
        for j in ti.static(range(4)):
            for k in ti.static(range(3)):
                w = w_center[i].x * w_side[j].y * w_center[k].z
                idx = (idx_center[i].x, idx_side[j].y, idx_center[k].z)
                grid_v_y[idx] += vp.y *w
                grid_w_y[idx] += w

    #for i, j, k in ti.ndrange(3, 3, 4):
    for i in ti.static(range(3)):
        for j in ti.static(range(3)):
            for k in ti.static(range(4)):
                w = w_center[i].x * w_center[j].y * w_side[k].z
                idx = (idx_center[i].x, idx_center[j].y, idx_side[k].z)
                grid_v_z[idx] += vp.z * w
                grid_w_z[idx] += w

代码文件:Reisen/3d_pic.py at main · mrzhuzhe/Reisen · GitHub

这个很简单,你只要在P2G和P2G这两个部分,根据particle位置计算grid index的时候考虑grid的平移就可以了。

其实[3,3,3]也是可以的吧。