作者：馬遙

最近胡淵鳴的文章“99行代碼的《冰雪奇緣》”掀起了一陣不大不小的熱度，而且還上了知乎熱榜。

傳送門：https://zhuanlan.zhihu.com/p/97700605

本來我只是簡單讀了一遍這篇文章，沒打算深究。但是由于這篇文章熱度很高，很多同學(xué)都問我這99行到底寫了些什么，該如何理解？特別是：如果要改動效果，該改哪些參數(shù)？

由于缺乏一點(diǎn)必要的基礎(chǔ)知識，新同學(xué)感覺無從下手。甚至有人因?yàn)閷?shí)在看不懂這99行代碼而打擊了自信心【笑】。

因此我研究了一下代碼，并寫出這篇技術(shù)分析的短文，同時(shí)加上詳細(xì)的代碼注釋。雖然仍然看不懂核心計(jì)算過程，但能讓新人至少理解這99行代碼的結(jié)構(gòu)，同時(shí)也能方便其它人參考，節(jié)約分析代碼的時(shí)間。

完整的代碼（含注釋）在文章末尾，可以對比查閱。

1、作為非專業(yè)領(lǐng)域的工程師，只能從程序設(shè)計(jì)角度觀其大略

首先說明，雖然很多人都是“職業(yè)程序員”，但是因?yàn)轭I(lǐng)域不同，實(shí)際工作內(nèi)容差別不小。比如做電商、網(wǎng)頁、電子游戲、手機(jī)應(yīng)用等等不同的領(lǐng)域，主要工作同樣是編程，工作的重點(diǎn)和知識體系都會有區(qū)別。

更不用說像原作者胡同學(xué)是從事算法相關(guān)的科研工作的，和一般的軟件工程師差別就更大了。

所以對這99行代碼不感興趣，或是因?yàn)橐恍╅T檻看不太明白也屬于正?，F(xiàn)象。比如說到“材料模擬”，就必然會有一堆Lame系數(shù)、μ、λ之類的東西，外行人當(dāng)然會直接蒙圈。

如果不熟悉python的numpy等科學(xué)運(yùn)算的庫，不熟悉矩陣運(yùn)算，那要想看明白就更需要花點(diǎn)功夫了。

好在一般的工程師也并不關(guān)心材料模擬算法的細(xì)節(jié)，我們關(guān)心的主要是：這個(gè)程序框架是怎樣的？如何修改？怎么用在實(shí)際場景里？ 等等這類更表面、更實(shí)際的問題。

2、基本數(shù)據(jù)類型說明

taichi庫最常用的數(shù)據(jù)類型是矩陣，ti.Vector、ti.Matrix、ti.var 這三個(gè)函數(shù)生成的都是矩陣，別被函數(shù)名騙了。具體說明：

1、ti.Vector，創(chuàng)建一個(gè)矩陣，矩陣的每個(gè)元素是向量。參數(shù)如下：

ti.Vector( 每個(gè)元素是幾維向量，dt=數(shù)據(jù)類型，shape=矩陣形狀）

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

矩陣的行數(shù)和列數(shù)由shape確定，如果shape是一個(gè)數(shù)字，就是單行矩陣。如果shape是一個(gè)元組(3,4)，就是3行4列矩陣。

例子：ti.Vector(2, dt=ti.f32, shape=1000)。它是1000個(gè)2維向量組成的1行矩陣，每個(gè)數(shù)字是float32。訪問最后一個(gè)元素寫v[999][1]即可

2、ti.Matrix，創(chuàng)建一個(gè)矩陣，矩陣的每個(gè)元素也是矩陣。參數(shù)如下：

ti.Matrix( 每個(gè)元素行數(shù)，每個(gè)元素列數(shù)，dt=數(shù)據(jù)類型，shape=矩陣形狀）

3、ti.Var，創(chuàng)建一個(gè)矩陣，每個(gè)元素是一個(gè)數(shù)值。參數(shù)如下：

ti.var(dt=數(shù)據(jù)類型, shape=矩陣形狀)

數(shù)據(jù)結(jié)構(gòu)決定處理方式，所以多花一點(diǎn)時(shí)間搞懂?dāng)?shù)據(jù)結(jié)構(gòu)，后面看算法就不容易暈了。比如代碼前面的矩陣定義：

x = ti.Vector(2, dt=ti.f32, shape=n_particles) # position位置, 每個(gè)粒子有一個(gè)位置

v = ti.Vector(2, dt=ti.f32, shape=n_particles) # velocity速度，每個(gè)粒子有一個(gè)速度

C = ti.Matrix(2, 2, dt=ti.f32, shape=n_particles) # affine速度場，每個(gè)粒子對應(yīng)一個(gè)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，這個(gè)例子里有3種材料，分別是0、1、2

Jp = ti.var(dt=ti.f32, shape=n_particles) # plastic deformation 塑性變形，不可恢復(fù)變形

grid_v = ti.Vector(2, dt=ti.f32, shape=(n_grid, n_grid)) # grid node momemtum/velocity 一個(gè)128x128矩陣，每個(gè)元素是一個(gè)Vector2。每個(gè)格子一個(gè)總體速度

grid_m = ti.var(dt=ti.f32, shape=(n_grid, n_grid)) # grid node mass格子質(zhì)量。128x128矩陣，每個(gè)元素是一個(gè)數(shù)字，每個(gè)格子一個(gè)質(zhì)量

?

關(guān)鍵是知道這個(gè)矩陣多大，是一個(gè)粒子對應(yīng)一個(gè)值，還是一個(gè)格子對應(yīng)一個(gè)值。

3、粒子particle，格子grid

熟悉游戲引擎優(yōu)化的同學(xué)應(yīng)該對“格子”這東西不陌生，2D游戲引擎有一種著名的優(yōu)化方式叫做“臟矩形”，和這次遇到的“格子grid”有相通之處。

99行程序的一開頭是這樣定義的：

# 粒子數(shù)量，網(wǎng)格數(shù)量=128*1

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

# 每個(gè)網(wǎng)格寬度ΔX=1/128，以及它的倒數(shù)inv_dx

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

?

（所有的坐標(biāo)、距離都是用0~1的小數(shù)表示）

可以看出，粒子總數(shù)很多，但是格子只有128*128個(gè)，示意圖大概是這樣：

示意圖，格子畫得很稀疏，實(shí)際上每個(gè)格子很小。上圖紅色圈是有粒子的格子，而藍(lán)色圈是沒有粒子的格子，藍(lán)色圈可以跳過計(jì)算，極大提高運(yùn)算效率。

原程序的思路是：某些算法與粒子緊密相關(guān)，每幀每個(gè)粒子都要計(jì)算；而另外一些屬性是總體屬性，只對每個(gè)包含粒子的格子計(jì)算一次。比如重力的影響就是整體屬性，要通過格子算。

而且，每個(gè)格子只對周圍相鄰的格子有影響，影響不會超過一個(gè)格子的距離。如果沒有這個(gè)假設(shè)，粒子的實(shí)時(shí)模擬就不可能做到了。

之后會看到，每一幀的計(jì)算分三段：

1. 粒子相關(guān)計(jì)算，粒子被歸入某一個(gè)格子（particle to grid）；

2. 格子之內(nèi)的計(jì)算，以及周邊一共9個(gè)格子互相的影響

3. 格子的速度、速度場等屬性，要再影響到每一個(gè)粒子（grid to particle）

4、taichi這個(gè)庫到底做了什么？

可以看出，所有粒子計(jì)算的方法，全都寫在了python代碼中，由此可見taichi這個(gè)庫并非一個(gè)即插即用的“物理模擬引擎”，而是一種用于科學(xué)計(jì)算的基礎(chǔ)設(shè)施~~

熟悉深度學(xué)習(xí)的朋友應(yīng)該對這種模式更熟悉一些，最常用的TensorFlow的底層也是在做同樣的事，實(shí)際上這些庫主要是解決了底層計(jì)算問題，要拿來直接用還得再做一層封裝。

★ 和深度學(xué)習(xí)的python代碼相同，程序并非是在執(zhí)行python腳本時(shí)進(jìn)行的計(jì)算，實(shí)際的底層運(yùn)算是被延后的。python腳本所做的事情都可以看成“準(zhǔn)備工作”。例如代碼中間這一行：

??? # 二次核函數(shù)

??? w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]

?

如果你想通過print 打印 w 的值只會無功而返，因?yàn)檫\(yùn)行完這句話，w的值還沒有開始計(jì)算呢~~這也為分析代碼帶來很多麻煩，因?yàn)椴荒芡ㄟ^打印log的方式觀察每一步計(jì)算的結(jié)果。（taichi自帶的print函數(shù)，我試過也沒成功，可能是方法不對）。

5、完整代碼+注釋

import taichi as ti

# 計(jì)算品質(zhì)，越大算得越準(zhǔn)確，但也越慢。

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

# 粒子數(shù)量=9000，網(wǎng)格數(shù)量=128

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

# 每個(gè)網(wǎng)格寬度ΔX=1/128，以及它的倒數(shù)inv_dx

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

# deltaTime，演算的時(shí)間間隔

dt = 1e-4 / quality

# 體積vol，密度 (rho就是ρ)

p_vol, p_rho = (dx * 0.5)**2, 1

# 質(zhì)量 = 體積 * 密度

p_mass = p_vol * p_rho

#以下是材料系數(shù)

# E=0.1e4=1000, 泊松分布系數(shù)nu=0.2

E, nu = 0.1e4, 0.2 # Young's modulus and Poisson's ratio

# Lame系數(shù)，定義材料性質(zhì)的，分別是μ和λ

mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1+nu) * (1 - 2 * nu)) # Lame parameters

?

?

# 小技巧：taichi的對象類型不易查看，可以調(diào)用矩陣對象的.to_numpy()方法，把它變成numpy對象再看

x = ti.Vector(2, dt=ti.f32, shape=n_particles) # position位置, 每個(gè)粒子有一個(gè)位置

v = ti.Vector(2, dt=ti.f32, shape=n_particles) # velocity速度，每個(gè)粒子有一個(gè)速度

C = ti.Matrix(2, 2, dt=ti.f32, shape=n_particles) # affine速度場，每個(gè)粒子對應(yīng)一個(gè)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，這個(gè)例子里有3種材料，分別是0、1、2

Jp = ti.var(dt=ti.f32, shape=n_particles) # plastic deformation 塑性變形，不可恢復(fù)變形

grid_v = ti.Vector(2, dt=ti.f32, shape=(n_grid, n_grid)) # grid node momemtum/velocity 一個(gè)128x128矩陣，每個(gè)元素是一個(gè)Vector2。每個(gè)格子一個(gè)總體速度

grid_m = ti.var(dt=ti.f32, shape=(n_grid, n_grid)) # grid node mass格子質(zhì)量。128x128矩陣，每個(gè)元素是一個(gè)數(shù)字，每個(gè)格子一個(gè)質(zhì)量

?

ti.cfg.arch = ti.cuda # Try to run on GPU，Nvidia顯卡的CUDA

?

?

# 下面的函數(shù)是每幀更新用的，先跳過這個(gè)函數(shù)看主程序初始化~~

# 留到最后看這個(gè)函數(shù)

@ti.kernel

def substep():

? # 每次都要將每個(gè)格子的總速度和總質(zhì)量歸0

? for i, j in ti.ndrange(n_grid, n_grid):

??? grid_v[i, j] = [0, 0]

??? grid_m[i, j] = 0

? # 更新每一個(gè)粒子，并讓它們屬于某一個(gè)格子(Particle to Grid)

? for p in range(n_particles): # Particle state update and scatter to grid (P2G)

??? # 用粒子坐標(biāo)x[p]，計(jì)算出它所屬的格子的左下角

??? base = (x[p] * inv_dx - 0.5).cast(int)

??? # fx是該粒子距離格子左下角的局部坐標(biāo)

??? fx = x[p] * inv_dx - base.cast(float)

?

??? # 二次核函數(shù)

??? # 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是一個(gè)延遲計(jì)算的表達(dá)式，無法直接看到它的值，后面很多參數(shù)都是類似的

??? # w是位置的函數(shù)。

??? # w是一個(gè)權(quán)重，會影響當(dāng)前格子與周圍格子的質(zhì)量、速度

??? # 猜測：一個(gè)格子邊緣的粒子，顯然對旁邊格子的影響更大

?

??? #每幀更新F，F(xiàn) = somefunc(F, C, dt)

??? F[p] = (ti.Matrix.identity(ti.f32, 2) + dt * C[p]) @ F[p] # deformation gradient update

?

??? # h "加硬"系數(shù)（雪積累在一起變硬），是根據(jù)塑性變形參數(shù)Jp算出來的

??? h = ti.exp(10 * (1.0 - Jp[p])) # Hardening coefficient: snow gets harder when compressed

??? # jelly 果凍的加硬系數(shù)固定0.3

??? if material[p] == 1: # jelly, make it softer

????? h = 0.3

??? # μ和λ的值根據(jù)h確定

??? mu, la = mu_0 * h, lambda_0 * h

??? # 液體的μ固定為0

??? if material[p] == 0: # liquid

????? mu = 0.0

??? # ---------------硬核計(jì)算開始，按理說不要隨意改動-----------------

??? U, sig, V = ti.svd(F[p])

??? J = 1.0

??? for d in ti.static(range(2)):

????? new_sig = sig[d, d]

????? if material[p] == 2:? # Snow

??????? new_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3)? # Plasticity 塑性

????? Jp[p] *= sig[d, d] / new_sig

????? sig[d, d] = new_sig

????? J *= new_sig

??? if material[p] == 0:? # Reset deformation gradient to avoid numerical instability

????? F[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 plasticity

??? stress = 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) * stress

??? affine = stress + p_mass * C[p]

??? # ===============硬核計(jì)算結(jié)束==================

?

??? # 遍歷相鄰8+1個(gè)格子，把粒子參數(shù)算到到格子上

??? for i, j in ti.static(ti.ndrange(3, 3)): # Loop over 3x3 grid node neighborhood

????? offset = ti.Vector([i, j])

????? dpos = (offset.cast(float) - fx) * dx

????? weight = w[i][0] * w[j][1]

????? # 每個(gè)粒子的速度、質(zhì)量疊加，得到當(dāng)前格子與周圍格子的速度、質(zhì)量

????? 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):

??? #如果這塊格子有質(zhì)量才需要計(jì)算

??? if grid_m[i, j] > 0: # No need for epsilon here

????? # 動量轉(zhuǎn)為速度。格子質(zhì)量越大，格子速度變得越小

????? grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] # Momentum to velocity

????? # 重力影響

????? grid_v[i, j][1] -= dt * 50 # gravity

????? # 碰到四周墻壁，格子速度強(qiáng)行置為0。數(shù)字3就是第幾個(gè)格子算墻壁的意思，可以改大試試

????? if i < 3 and grid_v[i, j][0] < 0:????????? grid_v[i, j][0] = 0 # Boundary conditions

????? if i > n_grid - 3 and grid_v[i, j][0] > 0: grid_v[i, j][0] = 0

????? if j < 3 and grid_v[i, j][1] < 0:????????? grid_v[i, j][1] = 0

????? if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0

?

?

? # 最后，把格子的計(jì)算還原到粒子上去。遍歷所有粒子

? 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)

??? # 同樣，要考慮該粒子周圍的幾個(gè)格子

??? for i, j in ti.static(ti.ndrange(3, 3)): # loop over 3x3 grid node neighborhood

????? dpos = ti.Vector([i, j]).cast(float) - fx

????? g_v = grid_v[base + ti.Vector([i, j])]

????? weight = w[i][0] * w[j][1]

????? # 新速度 += 權(quán)重 * 格子速度

????? new_v += weight * g_v

????? new_C += 4 * inv_dx * weight * ti.outer_product(g_v, dpos)

??? # 更新粒子的速度、C

??? v[p], C[p] = new_v, new_C

??? # 位置 += 速度 * ΔTime

??? x[p] += dt * v[p] # advection

?

# ~~~~初始化開始~~~~~

# 這里要結(jié)合運(yùn)行效果分析

import random

group_size = n_particles // 3?????? # 分三組，每組一個(gè)正方形

for i in range(n_particles):

? # 位置都是歸1化的，取值0~1

? x[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]???? # 初始速度0

? F[i] = [[1, 0], [0, 1]]?????? #變形梯度矩陣初始值

? Jp[i] = 1???????????????????? #塑料變形參數(shù)初始值

?

import numpy as np

# 初始化ti.GUI，顯示用

gui = ti.GUI("Taichi MLS-MPM-99", res=512, background_color=0x112F41)

# 運(yùn)行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)??? # 三種顏色，每組一個(gè)顏色

? gui.circles(x.to_numpy(), radius=1.5, color=colors[material.to_numpy()])??? # 根據(jù)位置畫小圓點(diǎn)，顏色3選1

? gui.show() # Change to gui.show(f'{frame:06d}.png') to write images to disk?? # 每幀畫一次

?

解析暫時(shí)先寫到這里。

本人畢竟不是專業(yè)研究物理模擬的學(xué)術(shù)界人士，深入閱讀代碼之后我自己也有一些疑問。我會根據(jù)讀者反饋繼續(xù)改進(jìn)這篇文章，能起到拋磚引玉的效果我就心滿意足了。

歡迎加入游戲開發(fā)群歡樂攪基：869551769?

有意向參與線下游戲開發(fā)學(xué)習(xí)的讀者可戳這里進(jìn)一步了解：http://www.levelpp.com/