taichi 冰雪奇缘学习

一、Taichi设计目标

性能

生产力

空间稀疏计算

可微编程 --> 深度学习

元编程

二、设计决策

计算与数据结构解耦合

自动特定的编译器优化

Megakernels

实现两个尺度自动微分

嵌入进python

三、99行代码分析

https://zhuanlan.zhihu.com/p/97700605

https://www.bilibili.com/read/cv4358206?spm_id_from=333.788.b_636f6d6d656e74.55

1、基本数据类型

taichi库最常用到的数据类型是矩阵，ti.var、ti.Vector、ti.Matrix这三个函数生成的其实都是矩阵，别被名字误导

1.1、ti.Vector

功能：创建一个矩阵，矩阵的每个元素是向量

ti.Vector(每个元素是几维向量，dt=数据类型，shape=矩阵形状)
dt：data_type的意思，dt=ti.f32指32位float
shape:确定矩阵的行数和列数，如果shape是一个数字，就是单行矩阵，如果shape是一个元祖（3,4），就是3行4列矩阵exp: ti.Vector(2,dt=ti.f32,shape=1000) 他是1000个2维向量组成的1行矩阵，每个数字是float32，访问最后一个元素写v[999][1]即可

1.2、ti.Matrix

功能：创建一个矩阵，矩阵的每个元素也是矩阵

ti.Matrix(每个元素行数，每个元素列数，dt=数据类型，shape=矩阵形状)

1.3、ti.var

功能：创建一个矩阵，每个元素是一个数值

ti.var(dt=数据类型，shape=矩阵形状)

1.4、数据结构分析

x=ti.Vector(2,dt=ti.f32,shape=n_particles) #position位置，每一个粒子有一个位置
v=ti.Vector(2,dt=ti.f32,shape=n_particles) #velocity速度，每一个粒子有一个速度
C=ti.Matrix(2,2,dt=ti.f32,shape=n_particles) #affine速度场，每一个粒子对应一个2X2矩阵
F=ti.Matrix(2,2,dt=ti.f32,shape=n_particles) #deformation gradient变形梯度矩阵
material=ti.var(dt=ti.i32,shape=n_particles) # material id 这个粒子里有三种材料，分别是0,1,2
Jp=ti.var(dt=ti.f32,shape=n_particles) # plastic deformation 塑形变形，不可恢复变形
grid_v=ti.Vector(2,dt=ti.f32,shape=(n_grid,n_grid)) #grid node
momemtum/velocity: 128X128矩阵,每一个元素是一个Vector2，每个格子一个总体速度
grid_m=ti.var(dt=ti.f32,shape=(n_grid,n_grid)) #grid node mass格子质量，128X128矩阵，每个元素是一个数字，每个格子一个质量

2、粒子particle，格子grid

# 99行程序开头定义
#粒子数量，网格数量=128*1
n_particles,n_grid=9000*quality**2,128*quality
#每个网格宽度dX=1/128,以及它的倒数inv_dx
dx,inv_dx=1/n_grid,float(n_grid)

所有的坐标、距离都是用0~1的小数表示

可以看出，粒子总数很多，但是格子只有128*128个，如下所示：

示意图中，格子画的很稀疏，实际上每个格子很小。红色圈是有粒子的格子，蓝色圈是没有粒子的格子，蓝色圈可以跳过计算，极大提高运算效率。

原程序的思路是：某些算法与粒子紧密相关，每帧每个粒子都要计算；而另外一些属性是总体属性，只对每个包含粒子的格子计算一次。比如重力的影响就是整体属性，要通过格子算。而且，每个格子只对周围相邻的格子有影响，影响不会超过一个格子的距离。如果没有这个假设，粒子的实时模拟就不可能做到了。

之后会看到，每一帧的计算分三段：

粒子相关计算，粒子被归入某一个格子（particle to grid）；
格子之内的计算，以及周边一共9个格子互相的影响
格子的速度、速度场等属性，要再影响到每一个粒子（grid to particle）

3、taichi这个库到底做了什么

可以看出，所有粒子计算的方法，全都写在了python代码中，由此可见taichi这个库并非一个即插即用的“物理模拟引擎”，而是一种用于科学计算的基础设施。

熟悉深度学习的朋友应该对这种模式更熟悉一些，最常用的TensorFlow的底层也是在做同样的事，实际上这些库主要是解决了底层计算问题，要拿来直接用还得再做一层封装。

和深度学习的python代码相同，程序并非是在执行python脚本时进行的计算，实际的底层运算是被延后的。python脚本所做的事情都可以看成“准备工作”。例如代码中间这一行：

# 二次核函数
w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]

如果你想通过print 打印 w 的值只会无功而返，因为运行完这句话，w的值还没有开始计算呢~~这也为分析代码带来很多麻烦，因为不能通过打印log的方式观察每一步计算的结果。

4、完整代码+注释

import taichi as ti# 计算品质，越大算得越准确，但也越慢。quality = 1 # Use a larger value for higher-res simulations# 粒子数量=9000，网格数量=128n_particles, n_grid = 9000 * quality ** 2, 128 * quality# 每个网格宽度ΔX=1/128，以及它的倒数inv_dxdx, inv_dx = 1 / n_grid, float(n_grid)# deltaTime，演算的时间间隔dt = 1e-4 / quality# 体积vol，密度 (rho就是ρ)p_vol, p_rho = (dx * 0.5)**2, 1# 质量 = 体积 * 密度p_mass = p_vol * p_rho#以下是材料系数# E=0.1e4=1000, 泊松分布系数nu=0.2E, nu = 0.1e4, 0.2 # Young's modulus and Poisson's ratio# Lame系数，定义材料性质的，分别是μ和λmu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1+nu) * (1 - 2 * nu)) # Lame parameters# 小技巧：taichi的对象类型不易查看，可以调用矩阵对象的.to_numpy()方法，把它变成numpy对象再看x = ti.Vector(2, dt=ti.f32, shape=n_particles) # position位置, 每个粒子有一个位置v = ti.Vector(2, dt=ti.f32, shape=n_particles) # velocity速度，每个粒子有一个速度C = ti.Matrix(2, 2, dt=ti.f32, shape=n_particles) # affine速度场，每个粒子对应一个2x2矩阵F = ti.Matrix(2, 2, dt=ti.f32, shape=n_particles) # deformation gradient变形梯度矩阵material = ti.var(dt=ti.i32, shape=n_particles) # material id，这个例子里有3种材料，分别是0、1、2Jp = ti.var(dt=ti.f32, shape=n_particles) # plastic deformation 塑性变形，不可恢复变形grid_v = ti.Vector(2, dt=ti.f32, shape=(n_grid, n_grid)) # grid node momemtum/velocity 一个128x128矩阵，每个元素是一个Vector2。每个格子一个总体速度grid_m = ti.var(dt=ti.f32, shape=(n_grid, n_grid)) # grid node mass格子质量。128x128矩阵，每个元素是一个数字，每个格子一个质量ti.cfg.arch = ti.cuda # Try to run on GPU，Nvidia显卡的CUDA# 下面的函数是每帧更新用的，先跳过这个函数看主程序初始化~~# 留到最后看这个函数@ti.kerneldef substep():# 每次都要将每个格子的总速度和总质量归0for i, j in ti.ndrange(n_grid, n_grid):grid_v[i, j] = [0, 0]grid_m[i, j] = 0# 更新每一个粒子，并让它们属于某一个格子(Particle to Grid)for p in range(n_particles): # Particle state update and scatter to grid (P2G)# 用粒子坐标x[p]，计算出它所属的格子的左下角base = (x[p] * inv_dx - 0.5).cast(int)# fx是该粒子距离格子左下角的局部坐标fx = x[p] * inv_dx - base.cast(float)# 二次核函数# Quadratic kernels  [http://mpm.graphics   Eqn. 123, with x=fx, fx-1,fx-2]w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]# w是一个延迟计算的表达式，无法直接看到它的值，后面很多参数都是类似的# w是位置的函数。# w是一个权重，会影响当前格子与周围格子的质量、速度# 猜测：一个格子边缘的粒子，显然对旁边格子的影响更大#每帧更新F，F = somefunc(F, C, dt)F[p] = (ti.Matrix.identity(ti.f32, 2) + dt * C[p]) @ F[p] # deformation gradient update# h "加硬"系数（雪积累在一起变硬），是根据塑性变形参数Jp算出来的h = ti.exp(10 * (1.0 - Jp[p])) # Hardening coefficient: snow gets harder when compressed# jelly 果冻的加硬系数固定0.3if material[p] == 1: # jelly, make it softerh = 0.3# μ和λ的值根据h确定mu, la = mu_0 * h, lambda_0 * h# 液体的μ固定为0if material[p] == 0: # liquidmu = 0.0# ---------------硬核计算开始，按理说不要随意改动-----------------U, sig, V = ti.svd(F[p])J = 1.0for d in ti.static(range(2)):new_sig = sig[d, d]if material[p] == 2:  # Snownew_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3)  # Plasticity 塑性Jp[p] *= sig[d, d] / new_sigsig[d, d] = new_sigJ *= new_sigif material[p] == 0:  # Reset deformation gradient to avoid numerical instabilityF[p] = ti.Matrix.identity(ti.f32, 2) * ti.sqrt(J)elif material[p] == 2:F[p] = U @ sig @ V.T() # Reconstruct elastic deformation gradient after plasticitystress = 2 * mu * (F[p] - U @ V.T()) @ F[p].T() + ti.Matrix.identity(ti.f32, 2) * la * J * (J - 1)stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stressaffine = stress + p_mass * C[p]# ===============硬核计算结束==================# 遍历相邻8+1个格子，把粒子参数算到到格子上for i, j in ti.static(ti.ndrange(3, 3)): # Loop over 3x3 grid node neighborhoodoffset = ti.Vector([i, j])dpos = (offset.cast(float) - fx) * dxweight = w[i][0] * w[j][1]# 每个粒子的速度、质量叠加，得到当前格子与周围格子的速度、质量grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)grid_m[base + offset] += weight * p_mass#遍历所有的格子for i, j in ti.ndrange(n_grid, n_grid):#如果这块格子有质量才需要计算if grid_m[i, j] > 0: # No need for epsilon here# 动量转为速度。格子质量越大，格子速度变得越小grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] # Momentum to velocity# 重力影响grid_v[i, j][1] -= dt * 50 # gravity# 碰到四周墙壁，格子速度强行置为0。数字3就是第几个格子算墙壁的意思，可以改大试试if i < 3 and grid_v[i, j][0] < 0:          grid_v[i, j][0] = 0 # Boundary conditionsif i > n_grid - 3 and grid_v[i, j][0] > 0: grid_v[i, j][0] = 0if j < 3 and grid_v[i, j][1] < 0:          grid_v[i, j][1] = 0if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0# 最后，把格子的计算还原到粒子上去。遍历所有粒子for p in range(n_particles): # grid to particle (G2P)base = (x[p] * inv_dx - 0.5).cast(int)fx = x[p] * inv_dx - base.cast(float)w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1.0), 0.5 * ti.sqr(fx - 0.5)]new_v = ti.Vector.zero(ti.f32, 2)new_C = ti.Matrix.zero(ti.f32, 2, 2)# 同样，要考虑该粒子周围的几个格子for i, j in ti.static(ti.ndrange(3, 3)): # loop over 3x3 grid node neighborhooddpos = ti.Vector([i, j]).cast(float) - fxg_v = grid_v[base + ti.Vector([i, j])]weight = w[i][0] * w[j][1]# 新速度 += 权重 * 格子速度new_v += weight * g_vnew_C += 4 * inv_dx * weight * ti.outer_product(g_v, dpos)# 更新粒子的速度、Cv[p], C[p] = new_v, new_C# 位置 += 速度 * ΔTimex[p] += dt * v[p] # advection# ~~~~初始化开始~~~~~# 这里要结合运行效果分析import randomgroup_size = n_particles // 3       # 分三组，每组一个正方形for i in range(n_particles):# 位置都是归1化的，取值0~1x[i] = [random.random() * 0.2 + 0.3 + 0.10 * (i // group_size), random.random() * 0.2 + 0.05 + 0.32 * (i // group_size)]material[i] = i // group_size # 0: fluid流体 1: jelly果冻 2: snow雪v[i] = [0, 0]     # 初始速度0F[i] = [[1, 0], [0, 1]]       #变形梯度矩阵初始值Jp[i] = 1                     #塑料变形参数初始值import numpy as np# 初始化ti.GUI，显示用gui = ti.GUI("Taichi MLS-MPM-99", res=512, background_color=0x112F41)# 运行20000帧for frame in range(20000):# s从0到18，每帧要算19次for s in range(int(2e-3 // dt)):substep()colors = np.array([0x068587, 0xED553B, 0xEEEEF0], dtype=np.uint32)    # 三种颜色，每组一个颜色gui.circles(x.to_numpy(), radius=1.5, color=colors[material.to_numpy()])    # 根据位置画小圆点，颜色3选1gui.show() # Change to gui.show(f'{frame:06d}.png') to write images to disk   # 每帧画一次
作者：皮皮关做游戏
https://www.bilibili.com/read/cv4358206?spm_id_from=333.788.b_636f6d6d656e74.55
出处： bilibili