普通软件工程师如何阅读“99行代码的《冰雪奇缘》”

2020-01-15 19:01 作者:皮皮关做游戏 0人读过 | 我要投稿

作者：马遥

最近胡渊鸣的文章“99行代码的《冰雪奇缘》”掀起了一阵不大不小的热度，而且还上了知乎热榜。

传送门：https://zhuanlan.zhihu.com/p/97700605

本来我只是简单读了一遍这篇文章，没打算深究。但是由于这篇文章热度很高，很多同学都问我这99行到底写了些什么，该如何理解？特别是：如果要改动效果，该改哪些参数？

由于缺乏一点必要的基础知识，新同学感觉无从下手。甚至有人因为实在看不懂这99行代码而打击了自信心【笑】。

因此我研究了一下代码，并写出这篇技术分析的短文，同时加上详细的代码注释。虽然仍然看不懂核心计算过程，但能让新人至少理解这99行代码的结构，同时也能方便其它人参考，节约分析代码的时间。

完整的代码（含注释）在文章末尾，可以对比查阅。

1、作为非专业领域的工程师，只能从程序设计角度观其大略

首先说明，虽然很多人都是“职业程序员”，但是因为领域不同，实际工作内容差别不小。比如做电商、网页、电子游戏、手机应用等等不同的领域，主要工作同样是编程，工作的重点和知识体系都会有区别。

更不用说像原作者胡同学是从事算法相关的科研工作的，和一般的软件工程师差别就更大了。

所以对这99行代码不感兴趣，或是因为一些门槛看不太明白也属于正常现象。比如说到“材料模拟”，就必然会有一堆Lame系数、μ、λ之类的东西，外行人当然会直接蒙圈。

如果不熟悉python的numpy等科学运算的库，不熟悉矩阵运算，那要想看明白就更需要花点功夫了。

好在一般的工程师也并不关心材料模拟算法的细节，我们关心的主要是：这个程序框架是怎样的？如何修改？怎么用在实际场景里？等等这类更表面、更实际的问题。

2、基本数据类型说明

taichi库最常用的数据类型是矩阵，ti.Vector、ti.Matrix、ti.var 这三个函数生成的都是矩阵，别被函数名骗了。具体说明：

1、ti.Vector，创建一个矩阵，矩阵的每个元素是向量。参数如下：

ti.Vector( 每个元素是几维向量，dt=数据类型，shape=矩阵形状）

其中：dt就是data_type的意思，dt=ti.f32 是指32位float，

矩阵的行数和列数由shape确定，如果shape是一个数字，就是单行矩阵。如果shape是一个元组(3,4)，就是3行4列矩阵。

例子：ti.Vector(2, dt=ti.f32, shape=1000)。它是1000个2维向量组成的1行矩阵，每个数字是float32。访问最后一个元素写v[999][1]即可

2、ti.Matrix，创建一个矩阵，矩阵的每个元素也是矩阵。参数如下：

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

3、ti.Var，创建一个矩阵，每个元素是一个数值。参数如下：

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

数据结构决定处理方式，所以多花一点时间搞懂数据结构，后面看算法就不容易晕了。比如代码前面的矩阵定义：

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、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矩阵，每个元素是一个数字，每个格子一个质量

关键是知道这个矩阵多大，是一个粒子对应一个值，还是一个格子对应一个值。

3、粒子particle，格子grid

熟悉游戏引擎优化的同学应该对“格子”这东西不陌生，2D游戏引擎有一种著名的优化方式叫做“脏矩形”，和这次遇到的“格子grid”有相通之处。

99行程序的一开头是这样定义的：

# 粒子数量，网格数量=128*1

n_particles, n_grid = 9000 * quality ** 2, 128 * quality

# 每个网格宽度ΔX=1/128，以及它的倒数inv_dx

dx, inv_dx = 1 / n_grid, float(n_grid)

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

可以看出，粒子总数很多，但是格子只有128*128个，示意图大概是这样：

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

原程序的思路是：某些算法与粒子紧密相关，每帧每个粒子都要计算；而另外一些属性是总体属性，只对每个包含粒子的格子计算一次。比如重力的影响就是整体属性，要通过格子算。

而且，每个格子只对周围相邻的格子有影响，影响不会超过一个格子的距离。如果没有这个假设，粒子的实时模拟就不可能做到了。

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

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

4、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的方式观察每一步计算的结果。（taichi自带的print函数，我试过也没成功，可能是方法不对）。

5、完整代码+注释

import taichi as ti

# 计算品质，越大算得越准确，但也越慢。

quality = 1 # Use a larger value for higher-res simulations

# 粒子数量=9000，网格数量=128

n_particles, n_grid = 9000 * quality ** 2, 128 * quality

# 每个网格宽度ΔX=1/128，以及它的倒数inv_dx

dx, 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.2

E, 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对象再看