Taichi Implicit MPM Solver

g1n0st · 2020 年12 月 7 日 08:01

最近看的几篇 mpm 相关的 paper 里面都是使用的隐式积分器，就趁着学习的机会看完了 Optimization integrator for large time steps，用 Taichi 实现了implicit MPM solver.

Demo

video1

video

实现

实现主要参考的还是16年蒋老师的 mpm-course 以及胡老师第二节课讲的隐式求解标准步骤，做完 particles to grid 以后构造：
h(\mathbf{v}^{n+1})=\mathbf{M} \mathbf{v}^{n+1} - \Delta t \mathbf{f}(\mathbf{x}^n+\Delta t \mathbf{v}^{n+1}) - \mathbf{M} \mathbf{v}^n=0
然后用 Newton-Raphson solver 迭代求解：
\mathbf{v}^{(i + 1)} = \mathbf{v}^{(i)}-(\frac{\partial h}{\partial v}(\mathbf{v}^{(i)} ))^{-1}h(\mathbf{v}^{(i)} )
隐式构造出 \frac{\partial h}{\partial v}(\mathbf{v}^{(i)} ) 矩阵后，用 conjugate gradient linear-solver 迭代求解每一步的线性方程组。最后做 grid to particles.

并实现了 Optimization integrator for large time steps 里的优化，
把求零点问题看成是求以下能量函数的最低点：
E(v_i)=\sum_{i}\frac{1}{2}m_{i}\vert\vert v_i - v_i^n\vert \vert^2+ \Phi(x_i^n+\Delta t v_i)
并用 line search 优化和修正最后的方向。
为耦合牛顿迭代和线性方程组求解，引入 relative tolerance 为 min(\frac{1}{2}, \sigma \sqrt{max(\vert\vert\nabla E\vert \vert), \tau})

用了维度无关编程和 ti.static() 以及稀疏数据结构，同时支持2维和3维的情形。
为了能在实时下运行，demo 中把 newton_tolerance 和 linear_solve_tolerance 分别设为 10^{-2} 和 10^{-6} .

感想

求 \frac{\partial f}{\partial v} 非常painful，自己尝试推了一天，最后参考了蒋老师的 ziran2020 库里面的求解式子。以及在解方程的时候，collision 情况的考虑需要自己仔细想一想。

Code

代码托管在 Github:

github.com

g1n0st/GAMES201/blob/master/extra/implicitMPM/implicit_mpm.py

import numpy as np
import taichi as ti

ti.init(arch = ti.cpu)

MATERIAL_PHASE_FLUID = 0
MATERIAL_PHASE_SOLID = 1

@ti.data_oriented
class IMPLICIT_MPM:
    def __init__(self,
    category = MATERIAL_PHASE_SOLID, dim = 2, dt = 0.00004, n_particles = 5000,  n_grid = 128, gravity = 9.8, gravity_dim = 1,
    p_rho = 1.0, E = 50000, nu = 0.4,
    newton_max_iterations = 5, newton_tolerance = 1e-5, line_search = True, line_search_max_iterations = 10, linear_solve_tolerance_scale = 1,
    linear_solve_max_iterations = 1000, linear_solve_tolerance = 1e-30,
    cfl = 0.4, ppc = 16, 
    real = float,
    debug_mode = True,
    ignore_collision = False
    ):

This file has been truncated. show original

Future work

修复一些已知的 Bugs，改进数据布局和一些 quick-and-dirty 的实现。
由于要考虑 \Psi ,\frac{\partial \Psi}{\partial \mathbf{F}} ,\frac{\partial^2 \Psi}{\partial \mathbf{F} \partial \mathbf{F}} 等众多量的求解，目前只实现了一个 fixed corotated 模型，以后会抽象成一个专门的 stress 类。
用 possion disk sampling 实现 mesh to particles.
把 implicit solver 用于 IQMPM，anisoMPM 等。

References

[1] Gonzalez, O. and Stuart, A. (2008). A first course in continuum mechanics. Cambridge University Press.
[2] Gast, T., Schroeder, C., Stomakhin, A., Jiang, C., and Teran, J. (2015). Optimization integrator for large time steps. IEEE Trans Vis Comp Graph, 21(10):1103–1115.
[3] C. Jiang, C. Schroeder, J. Teran, A. Stomakhin, and A. Selle. 2016. The material point method for simulating continuum materials. In SIGGRAPH Course. 24:1–24:52.
[4] C. Jiang, C. Schroeder, A. Selle, J. Teran, and A. Stomakhin. 2015. The affine particle-in-cell method. ACM Trans Graph 34, 4 (2015), 51:1–51:10.

g1n0st · 2020 年12 月 9 日 13:01

目前line search的部分还存在bug，虽然residual减少了但是total energy反而增加了……还在debug中

mapengfei · 2021 年7 月 27 日 07:59

同学你好，我想问你一个牛顿迭代的问题。我看了你列的参考文献中第三本参考书，找到了公式（1）。

我对于牛顿迭代的步骤还不太清楚。假设在第k步牛顿迭代，我们能求出网格节点 i 上的速度为 u^{(k)}_{i} ，我们用这个速度去更新粒子的形变梯度，计算出 \frac{\partial^{2} \Psi}{\partial \mathrm{F} \partial \mathrm{F}}\left(\hat{\mathrm{F}}_{\mathrm{p}}(\hat{\boldsymbol{x}})\right) ，还要用这个速度更新粒子位移，才能进行P2G插值，最后求出 \delta \mathbf{f}_i ，这样对吗？

\begin{array}{c} -\delta \mathbf{f}_{\mathfrak{i}}=\sum_{j} \frac{\partial^{2} e}{\partial \hat{x}_{i} \partial \hat{x}_{j}}(\hat{x}) \delta u_{j}=\sum_{p} V_{p}^{0} A_{p}\left(F_{p}^{n}\right)^{\top} \nabla w_{i p}^{n} \\ A_{p}=\frac{\partial^{2} \Psi}{\partial F \partial F}\left(\hat{\mathbf{F}}_{p}(\hat{x})\right):\left(\sum_{j} \delta u_{j}\left(\nabla w_{j p}^{n}\right)^{\top} F_{p}^{n}\right) \end{array}\tag{1}