HW 0.5: 后向欧拉法求解弹簧质点系统

弹簧刚度超过 100 万,我们的弹簧还不爆炸。显示求解早炸了。
应该是成了

ms_implicit.py

import taichi as ti

# ti.init(debug=True, arch=ti.cpu)
ti.init(arch=ti.gpu)

## 系统参数 ===============================
MAX_NUM_PARTICLES   = 256   # 最多质点数
PARTICLE_MASS       = 1     # 质点质量
BOTTOM_Y            = 0.05  # 地面位置
CONNECTION_RADIUS   = 0.15  # 质点自动连接半径
GRAVITY = ti.Vector([0, -9.8], dt=ti.f32) # 重力场
dt = 1e-3                   # 时间步长

num_particles       = ti.var(ti.i32, shape=())  # 现有质点数
spring_stiffness    = ti.var(ti.f32, shape=())  # 弹簧刚度
paused              = ti.var(ti.i32, shape=())  # 暂停
damping             = ti.var(ti.f32, shape=())  # 阻尼
# rest_length[i, j] = 0 则 i,j 未连接
rest_length = ti.var(ti.f32, shape=(MAX_NUM_PARTICLES, MAX_NUM_PARTICLES))

x = ti.Vector(2, dt=ti.f32, shape=MAX_NUM_PARTICLES) # 位置
v = ti.Vector(2, dt=ti.f32, shape=MAX_NUM_PARTICLES) # 速度

M = ti.Matrix(2, 2, dt=ti.f32, shape=(MAX_NUM_PARTICLES, MAX_NUM_PARTICLES))
J = ti.Matrix(2, 2, dt=ti.f32, shape=(MAX_NUM_PARTICLES, MAX_NUM_PARTICLES))
A = ti.Matrix(2, 2, dt=ti.f32, shape=(MAX_NUM_PARTICLES, MAX_NUM_PARTICLES))
F = ti.Vector(2, dt=ti.f32, shape=MAX_NUM_PARTICLES)
b = ti.Vector(2, dt=ti.f32, shape=MAX_NUM_PARTICLES)

spring_stiffness[None] = 10000
damping[None] = 20  # 恒定阻尼


## 速度求解器 =====================================
@ti.kernel
def symplectic_euler():
    "半隐式欧拉法"
    n = num_particles[None]
    k = spring_stiffness[None]
    m, g  = PARTICLE_MASS, GRAVITY

    for i in range(n):
        v[i] *= ti.exp(-dt * damping[None]) # 计算阻尼
        f = m * g
        for j in range(n):
            l_ij = rest_length[i, j]
            if l_ij != 0:  # 两质点间有弹簧连接
                x_ij = x[i] - x[j]
                f += -k * (x_ij.norm() - l_ij) * x_ij.normalized()
        
        # `dv = dt * a = dt * (f / m)`
        v[i] += dt * (f / m)


## 隐式欧拉法 --------------
@ti.kernel
def init_M():
    """初始化质量矩阵 M
    """
    m = ti.Matrix([
        [PARTICLE_MASS, 0],
        [0, PARTICLE_MASS]
    ])
    for i in range(num_particles[None]):
        M[i, i] = m


@ti.kernel
def update_J():
    """更新 jacobi 矩阵
    - [Miles Macklin](https://blog.mmacklin.com/2012/05/04/implicitsprings/)
    """
    I = ti.Matrix([
        [1.0, 0.0],
        [0.0, 1.0]
    ])
    k = spring_stiffness[None]
    for i, d in J:  # 遍历 J
        # i 为观察的质点
        # d 为求导数的方向 x_d
        for j in range(num_particles[None]):  # 遍历所有点
            l_ij = rest_length[i, j] # 对于弹簧 i <-- j
            if (l_ij != 0) and (d == i or d == j):
                # i,j 间有弹簧链接 && 求导方向 d 为 i 或 j
                x_ij = x[i] - x[j]
                X_ij_bar = x_ij / x_ij.norm()
                mat = X_ij_bar.outer_product(X_ij_bar)
                J[i, d] += -k * (I - l_ij/x_ij.norm() * (I - mat))
                if d==i:
                    J[i, d] *=  1.0
                else: # d==j
                    J[i, d] *= -1.0


@ti.kernel
def update_A(beta: ti.f32):
    """更新 A
        A = M - dt^2 * J(t)
    """
    # beta = 0.5
    for i, j in A:
        A[i, j] = M[i, j] - beta * dt**2 * J[i, j]


@ti.kernel
def update_F():
    """计算 x 的受到的合力
    """
    k = spring_stiffness[None]
    m, g = PARTICLE_MASS, GRAVITY

    for i in range(num_particles[None]):
        F[i] = m * g

    for i, j in rest_length:
        l_ij = rest_length[i, j]
        if l_ij != 0:
            x_ij = x[i] - x[j]
            F[i] += -k * (x_ij.norm() - l_ij) * x_ij.normalized()
        else:  # l_ij == 0
            pass # i,j 间无弹簧


@ti.kernel
def update_b():
    """更新 b

        b = M*v_star + dt * F(t)
        v_star = v(t) + dt * a(t)

    `a(t)` 其他力导致的加速度。如:damping、相对速度
    """
    for i in range(num_particles[None]):
        v_star = v[i] * ti.exp(-dt * damping[None])
        b[i] = A[i, i] @ v_star + dt * F[i]


@ti.kernel
def jacobi():
    """Jacobi 迭代

        A = M - dt^2 * J(t)
        b = M * v(t) + dt * F(t)
        A * v(t+1) = b
    """
    n = num_particles[None]
    for i in range(n):
        for j in range(n):
            if i != j:
                b[i] -= A[i, j] @ v[j]

        v[i] = A[i, i].inverse() @ b[i]


def implicit_euler(beta=0.5):
    """隐式欧拉法 + Jacobi 迭代

        A = M - beta * dt^2 * J(t)
        b = M*v(t) + dt * F(t)
        A * v(t+1) = b

    ### beta
        = 0.0: forward/semi-implicit Euler (explicit)
        = 0.5: middle-point (implicit)
        = 1.0: backward Euler (implicit)

    ### step
    0. 初始化 M
    1. 更新 J(t)
    2. 更新 A
    3. 更新 F(t)
    4. 更新 b
    5. 求解 v(t+1)
    """
    init_M()
    update_J()
    update_A(beta)
    update_F()
    update_b()
    jacobi()


## 通用步骤 ====================================
@ti.kernel
def collide():
    """与地面碰撞"""
    for i in range(num_particles[None]):
        if x[i][1] < BOTTOM_Y:
            x[i][1] = BOTTOM_Y
            v[i][1] = 0


@ti.kernel
def update_position():
    """更新位置
        x(t+1) = x(t) + dt * v(t+1)
    """
    for i in range(num_particles[None]):
        x[i] += v[i] * dt


def substep():
    "一个时间步长"
    # symplectic_euler() # 半隐式欧拉法
    implicit_euler(1.0)    # 隐式欧拉法/后向欧拉法
    collide()
    update_position()


@ti.kernel
def add_particle(pos_x: ti.f32, pos_y: ti.f32):
    "添加新质点"
    new_particle_id = num_particles[None]
    x[new_particle_id] = [pos_x, pos_y]
    v[new_particle_id] = [0, 0]
    num_particles[None] += 1
    
    # Connect with existing particles
    for i in range(new_particle_id):
        dist = (x[new_particle_id] - x[i]).norm()
        if dist < CONNECTION_RADIUS: # 与指定半径内的点建立连接
            rest_length[i, new_particle_id] = 0.1
            rest_length[new_particle_id, i] = 0.1


## GUI =====================================
gui = ti.GUI('Mass Spring System', res=(512, 512), background_color=0xdddddd)

add_particle(0.3, 0.3)
add_particle(0.3, 0.4)
add_particle(0.4, 0.4)

while True:
    for e in gui.get_events(ti.GUI.PRESS):
        if e.key in [ti.GUI.ESCAPE]:
            exit()
        elif e.key == gui.SPACE:
            paused[None] = not paused[None]
        elif e.key == ti.GUI.LMB:
            add_particle(e.pos[0], e.pos[1])
        elif e.key == 'c':
            num_particles[None] = 0
            rest_length.fill(0)
        elif e.key == 's':
            if gui.is_pressed('Shift'):
                spring_stiffness[None] /= 1.1
            else:
                spring_stiffness[None] *= 1.1
        elif e.key == 'd':
            if gui.is_pressed('Shift'):
                damping[None] /= 1.1
            else:
                damping[None] *= 1.1

    if not paused[None]:
        for _ in range(10): # 10 小步为一帧
            substep()
    
    ## 绘制
    # 地板
    gui.line(begin=(0.0, BOTTOM_Y), end=(1.0, BOTTOM_Y), color=0x0, radius=1)

    # ms 系统
    X = x.to_numpy()
    gui.circles(X[:num_particles[None]], color=0xffaa77, radius=5)
    for i in range(num_particles[None]):
        for j in range(i + 1, num_particles[None]):
            if rest_length[i, j] != 0:
                gui.line(begin=X[i], end=X[j], radius=2, color=0x445566)
    
    ## 显示系统参数
    gui.text(content=f'C: clear all; Space: pause', pos=(0, 0.95), color=0x0)
    gui.text(content=f'S: Spring stiffness {spring_stiffness[None]:.1f}', pos=(0, 0.9), color=0x0)
    gui.text(content=f'D: damping {damping[None]:.2f}', pos=(0, 0.85), color=0x0)
    gui.show()

Jacobi 矩阵

根据 Miles Macklin 给出的公式稍微化简了下

image

equ LaTeX
$$
\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_i}=-k\left[ \mathbf{I}-\frac{l_{ij}}{\lVert \mathbf{x}_{ij} \rVert}\left( \mathbf{I}+\widehat{\mathbf{x}}_{ij}\cdot \widehat{\mathbf{x}}_{ij}^{\text{T}} \right) \right] 
\\
\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_j}=-\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_i}
\\
\mathbf{J}_{ij}=\sum_{\text{k}}^{\text{n}}{\frac{\partial \mathbf{f}_{ik}}{\partial \mathbf{x}_j}}
$$

跑起来感觉没问题
(看上去是写对了,jacobi 矩阵的推导先咕咕咕


我听说有个方法叫 Jacobian-free Newton Krylov,那可太棒了。坐等 lec 04

求导是不可能求的,只能看别人的结果推一推过程这样。

11 Likes

你好,上课的时候有个疑问想向你请教一下,为什么要先与地面碰撞然后再更新位置? 反过来,把与地面碰撞放在更新位置之后,不就能保证新的位置肯定不在bottom_y以下了吗?(我看在collide函数里面,会判断y是否小于bottom_y)

你这么一说,我感觉是有点问题。

貌似讲反了。

    """与地面碰撞"""
    for i in range(num_particles[None]):
        if x[i][1] < BOTTOM_Y:
            x[i][1] = BOTTOM_Y
            v[i][1] = 0

    """更新位置"""
    for i in range(num_particles[None]):
        x[i] += v[i] * dt

目前是先碰撞,在计算位置。
如果例子恰好在边界上(x[i][1] == BOTTOM_Y),而速度又为负数。则不会触发碰撞。
等到计算位置时,就跑到地底下去了。当然下一次循环就又被拉回来了。
但是还是有机会被看见。

反过来先计算位置。这时可能跑到地面以下,然后再碰撞,把它纠正过来。
看起来很合理,每一次迭代都不会导致粒子跑到地面以下。

不太理解A和b的表达式,跟课上给出的公式不一样

简单地说我方程两边都乘上了 M,矩阵相乘比求逆好一点,虽然这个就是个稀疏的对角矩阵。

1 Like

[Taichi] version 0.6.5, supported archs: [cpu, cuda, opengl], commit c3b7ce52, python 3.7.4
运行这个代码报错:
106 x_ij = x[i] - x[j]
107 X_ij_bar = x_ij / x_ij.norm()
–> 108 mat = X_ij_bar.outer_product(X_ij_bar)
109 J[i, d] += -k * (I - l_ij/x_ij.norm() * (I - mat))
110 if d==i:

TypeError: outer_product() missing 1 required positional argument: ‘b’
这个是以为api更改吗,还是我的版本不够新?

版本有点老了,更新到 0.6.9 吧

我有疑问,jacobi() 函数应该是迭代计算的吧?运行一次肯定不够的,至少运行几十次才能出正确结果吧?
像这样:

    update_F()
    update_b()
    for iii in range(50):
        jacobi()

另外,楼主的jacobi 函数似乎有点问题,缺少中间变量,多次运行后就飞了。我按照jacobi的范例,加上了中间变量new_v和r,结果似乎是正确了。

但是现在问题是,改为迭代运算后,运行极为缓慢,不知如何解决。
总之感谢楼主的偏导方程……这个问题的数学难度真的很大

1 Like

更新b的时候b[i] = A[i, i] @ v_star + dt * F[i]是不是错了?应该是b[i] = M[i, i] @ v_star + dt * F[i]?

计算Jacobian修改一下,以及按楼上说的rhs改成乘M,再把后端换成cpu,基本正常了

@ti.kernel
def update_J():
    """更新 jacobi 矩阵
    - [Miles Macklin](https://blog.mmacklin.com/2012/05/04/implicitsprings/)
    """
    I = ti.Matrix([
        [1.0, 0.0],
        [0.0, 1.0]
    ])
    k = spring_stiffness[None]
    for i, d in J:  # 遍历 J
        # i 为观察的质点
        # d 为求导数的方向 x_d
        J[i, d] *= 0.0
        for j in range(num_particles[None]):  # 遍历所有点
            l_ij = rest_length[i, j] # 对于弹簧 i <-- j
            if (l_ij != 0) and (d == i or d == j):
                # i,j 间有弹簧链接 && 求导方向 d 为 i 或 j
                x_ij = x[i] - x[j]
                X_ij_bar = x_ij / x_ij.norm()
                mat = X_ij_bar.outer_product(X_ij_bar)
                if d==i:
                    J[i, d] += -k * (I - l_ij/x_ij.norm() * (I - mat))
                else: # d==j
                    J[i, d] += k * (I - l_ij/x_ij.norm() * (I - mat))

周末得抽时间来研究一下哈哈

题主有一处笔误,latex简化公式的结果应该是,贴出来的公式最后一个减号写加号了:

\frac{\partial \mathbf{f}_{ij}}{\partial \mathbf{x}_i}=-k\left[ \mathbf{I}-\frac{l_{ij}}{\lVert \mathbf{x}_{ij} \rVert}\left( \mathbf{I}-\widehat{\mathbf{x}}_{ij}\cdot \widehat{\mathbf{x}}_{ij}^{\text{T}} \right) \right]

按照链接的公式复现了一遍,100w也是正常的