行列

内積

二次元ベクトルv・w の内積は、各成分の掛け算の和。
やっていることは、wの射影をvの長さでかけている。
これは逆にvの射影をwの長さでかけても同じ。

また、各成分の掛け算の和がどうして射影と関係あるのか？また次回

three.js

発散(divergence)

昨日、advect+forceで速度場にマウスの力を与えた。
可視化したら、マウスの範囲だけで少し流れるようなものが見えた。

マウスの力をあたえた結果、物理的に不正な値が生まれる。
マウスの中心には、「湧き出し」が発生。
水は圧縮できないので、これは正しくない。

divergence.fragでは、ただ湧き出しの量を求める。後の圧力(computePressure.frag)で発散を0にするため頑張る。

空間微分で解いていく。
空間微分とは有限差分法的に言うと隣同士のセルの差のことです。自分のセルの値は使わずに両隣の値の差分をとる中心差分の方がより精密になるみたいです。移流項以外をこの方法で解いていきます。
ナビエ・ストークス方程式に何度も出てくる「ナプラ」「▽」は空間微分のこと。
発散: 速度場の中心差分で ∂u/∂x + ∂v/∂y

発散はping-pong が不要なため、レンダーターゲットで対応。

uniform sampler2D uVelocity
uniform vec2 uPx // 隣の速度を読むのに使う
uniform float uDt // 1秒あたりの発散を求めるのに使う

void main() {
 float x0 = texture2D(uVelocity, vUv - vec2(uPx.x, 0.0)).x;
 float x1 = texture2D(uVelocity, vUv + vec2(uPx.x, 0.0)).x;
 float y0 = texture2D(uVelocity, vUv - vec2(0.0, uPx.y)).y;
 float y1 = texture2D(uVelocity, vUv + vec2(0.0, uPx.y)).y;

 float div = (x1 - x0) / 2.0 + (y1 - y0) / 2.0;

 gl_FragColor = vec4(div / max(uDt, 1e-6), 0.0, 0.0 , 1.0); //1秒あたりの発散
}

圧力

ヤコブ反復でpressureを求める。

まず、左右上下のpressureを読む。
同じ位置のdivergenceを読む。
newPressure=平均圧力 - divergence
ポアソン方程式の離散解法（Jacobi）1回分

これを複数回繰り返して pressure を収束させる。
具体的にはupdate()関数でループ。GPUComputationRendererを使っているのでcompute()する。

経過

computeVelocity 速場に移流＋マウスの力をプラス。
divergence 湧き出しを計算
computePressure 湧き出し/吸い出しを0にする圧力を求める。(x24)
圧力で速度を修正して、正しい速度場を求める。(発散ゼロに)

ここまで、1でもとめたuVelocityはまだ不正。

Projection

現在の速度を用意。

計算した圧力から、圧力の勾配を計算する。また空間微分（中間差分）を使う。
圧力勾配: 圧力場の中心差分で ∇p = (∂p/∂x, ∂p/∂y)

vec2 gradP = vec2(p0 - p1, p2 - p3) * 0.5;
v = v - gradP * uDt;

ヘルムホルツ分解の定理
任意のベクトル場vは以下のように分解できる。
v = v_div_free + grad(p)
v_div_free 発散ゼロの成分（今求めたいもの！）
つまりv_div_free = v - grad(p)

ここの*uDtは、divergence.fragで div / uDtをしているから。
ちなみに、両方のuDtを消すと、上手く機能しているように見えた。

これで、発散がゼロの正しい速度場が得られた。

gl_FragCoord = vec4(v, 0.0, 1.0);

また、この時点で流体はきれいに可視化可能。ここから色付。

dye

velocityとほとんど同じ。色が早く消えるほうが視覚的にきれいなため、velocityの減衰より低い(強い)減衰を渡す。

ラスト。出力。

内省

11月にやったときは1週間くらいじっと見て挫折したが、今GPUComputationRendererでリベンジしたら2日でラストまで行けて、成長を感じた。これは嬉しいなあ。

でも早くこれを理解して、ポートフォリオ作って、就職口探す。

But I need to understand this, make portfolio, and I start looking for a job.

BLOG

3/12

3/11

3/10

3/5

3/4

3/3

3/2

3/1

2/28

2/27