JosJos的平流項

Jos的程式碼裡，平流項看似簡單

GPU Gems 38
半拉格朗日法

幾何上算容易理解

但

為什麼不是直接使用 平流項(v·∇)v 呢？

I know nothing

お元気ですか？

下面就讓天添老師娓娓道來 

(轉身換手爬出場外) 

太極圖形課S1第11講：流體仿真 02
直接用平流項的有限差分
步長比較大，誤差的累積
會導致數值不穩定？

太極圖形課S1第11講：流體仿真 02
重整一下半拉格朗日法
現出原形
原來就是步長更短的有限差分

Fluid這條路上，爬著也算爬完了嗯 🤠🙃

回顧

這是NS方程

NS方程

這是Jos的Operator splitting

Operator splitting
招式表出處

是不是急著放學？ 🤠

老師再見 同學再見 😋🕹️

NS的平流項 (v·∇)v

(v·∇)              (∇·v)

欸不是嗎 🙃

太極圖形課S1第11講：流體仿真 02

(v·∇)v其實就是v·∇v

原來是括號之亂 😋🤠 奇怪的知識又增加了

JosJos的奇妙冒險

看不懂的算法 是否來自太空

別人寫的代碼 看著像密碼 🎶

    ── Booming Tech。擎空

事情是這樣的

幾年前我買了一本

Jos Stam的The Art of Fluid Animation中譯本

隨書附上的 程式碼

陸陸續續讀了3遍

有2個地方就是看不懂

和書的內容對不上，分別是

1. 壓力項：project函數
2. 擴散項：diffuse函數 

直到昨天看了

• GAMES103第11講
• GAMES201第4講

終於讓我找到關鍵的One Piece 🧐

拉普拉斯算子的離散形式

GAMES103第11講

projection函數在做什麼？

GAMES201第4講
上面先假設存在1個不可壓縮速度場u*
推導出了泊松方程

GAMES201第4講
這是上面泊松方程的離散形式
解泊松方程可以得出壓力場P

藍色的部分我先不管
Jos的P和u都已經攤平成1D陣列
要找出矩陣A滿足粉紅色的算式相對簡單
A x = b

Jos 在這裡做了一些優化
可以忽略Δt，因為那只是1個scale
當有一個線性方程組為
A x = (1/Δt) b
忽略上式右邊的 (1/Δt)後，會有
A (Δt x) = b
忽略(1/Δt)得到的結果相當於預乘了Δt

GAMES201第4講
有了P，就能算出不可壓縮的速度場u*

上面2個步驟對映的程式碼

1. 解泊松方程可以得出壓力場P
2. 有了P，就能算出不可壓縮的速度場u*

解泊松方程迭代的過程中

會使用到已更新的鄰居值

乍看之下會懷疑是不是程式碼寫錯了？

但它其實是在做 Gauss-Seidel迭代

Jacobi迭代

就比較接近一般Programmer的想法了

diffuse函數在做什麼？

diffuse 和 viscosity 有關

關於擴散項，一開始我們有

(ρ是流體密度，v = 黏滯系數/ρ)

 a = v▽·▽u  

u 是當前速度

u' 是下一步的速度

u'' 是下下一步的速度

用u時，當前會影響下一步

Δu = v▽·▽u 

u' - u = v▽·▽u 

用u'時，下一步會影響下下一步

Δu' = v▽·▽u' 

u'' - u' = v▽·▽u' 

上式移項後有

u'' = u' + v▽·▽u' 

既然往前+v▽·▽u' 會得到 u''

反過來 -v▽·▽u' 就會得到u吧？

u = u' - v▽·▽u' 

u' - u = v▽·▽u' 

Δu = v▽·▽u' 

這樣就和Jos的程式碼對上了

Δu = v▽·▽u'

對映的程式碼

結語

 

此碼只應天上有？ 🤠🤔 人間能得幾回聞 

補充資料

NS方程的幾何意義

From wiki

NS方程超清解說
太極圖形课S1第10講：流體仿真 01