Houdiniでぽあそん

上記のグラフの接続性は、位置(3,3) におけるグラフの接続性ですが、他のすべての点についても同様に考えられます。

\( L_{32} = \begin{matrix} [0]_{51} & - & [0]_{52} & - & [0]_{53} & - & [0]_{54} & - & [0]_{55} \\ | & &   | & &   | & & | & & | \\ [0]_{41} & - & [1]_{42} & - & [0]_{43} & - & [0]_{44} & - & [0]_{45} \\ | & &   | & &   | & & | & & | \\ [1]_{31} & - & [-4]_{32} & - & [1]_{33} & - & [0]_{34} & - & [0]_{35} \\ | & &   | & &   | & & | & & | \\ [0]_{21} & - & [1]_{22} & - & [0]_{23} & - & [0]_{24} & - & [0]_{25} \\ | & &   | & &   | & & | & & | \\ [0]_{11} & - & [0]_{12} & - & [0]_{13} & - & [0]_{14} & - & [0]_{15} \\ \end{matrix} \)

境界部分の点では境界条件を考える必要が出てきます。

これもポアソン方程式を解くうえで重要な要素となりますが、今回説明は割愛します。（値0のデイリクレ条件（固定境界）で実装していきます。）

\( L_{11} = \begin{matrix} [0]_{51} & - & [0]_{52} & - & [0]_{53} & - & [0]_{54} & - & [0]_{55} \\ | & &   | & &   | & & | & & | \\ [0]_{41} & - & [0]_{42} & - & [0]_{43} & - & [0]_{44} & - & [0]_{45} \\ | & &   | & &   | & & | & & | \\ [0]_{31} & - & [0]_{32} & - & [0]_{33} & - & [0]_{34} & - & [0]_{35} \\ | & &   | & &   | & & | & & | \\ [1]_{21} & - & [0]_{22} & - & [0]_{23} & - & [0]_{24} & - & [0]_{25} \\ | & &   | & &   | & & | & & | \\ [-4]_{11} & - & [1]_{12} & - & [0]_{13} & - & [0]_{14} & - & [0]_{15} \\ \end{matrix} \)

各点ごとのグラフの接続性が行（対称行列になるので列でもOK）となるように行列を構築することで、ラプラシアン行列を作成できます。

\( A = \begin{bmatrix} -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & ← L_{11} \\ \vdots & \ddots & & & & & & & & & & & & & & & & & & & & & & & \vdots & ← L_{..} \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & ← L_{32} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & ← L_{33} \\ \vdots & & & & & & & & & & & & & & & & & & & & & & & \ddots & \vdots & ← L_{..} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 &     -4 & ← L_{55} \\ \end{bmatrix} \)

上記のラプラシアン行列 A と解の値である f とでポアソン方程式を解くことができます。

\( A\phi = f \)

Houdiniで実装

上記をHoudiniで実装していきます。

1. 2次元の格子の代わりとして2Dボリュームを作成します。
   

   

2. Aφ = f の f の部分を適当に作成します。（電荷密度などに相当するもの）
   

   

3. Python SOPで Aφ = f の A の行列を作成し、np.linalg.solveを使い Aφ = f の φ について解きます。
   

   

   直接法 ソースを展開 ソースを折りたたむ

   ---

   node = hou.pwd()
   geo = node.geometry()
    
   # Add code to modify contents of geo.
   # Use drop down menu to select examples.
   import numpy as np
    
   volume = geo.prim(0)
   res = volume.resolution()
   f   = volume.allVoxels()
    
   A = np.zeros((res[0]*res[1], res[0]*res[1]))
   f = np.array(f)
   for y in range(res[1]):
       for x in range(res[0]):
           idx = x + y*res[0]
           A[idx,idx] = -4
           if x-1 > -1:
               idx1 = (x-1) + y*res[0]
               A[idx,idx1] = 1
           if x+1 < res[0]:
               idx2 = (x+1) + y*res[0]
               A[idx,idx2] = 1
           if y-1 > -1:
               idx3 = x + (y-1)*res[0]
               A[idx,idx3] = 1
           if y+1 < res[1]:
               idx4 = x + (y+1)*res[0]
               A[idx,idx4] = 1
          
   phi = np.linalg.solve(A, f)
   volume.setAllVoxels(phi)
   volume.setAttribValue("name", "phi")
   

   ---

4. φ について可視化し結果を確認します。
   

   

5. φ のラプラシアンに関しても可視化し、ほぼ f と一致していることが確認できます。
   よって、 Aφ = f  が問題なく解けていることが確認できます。
   
   

   左：φ のラプラシアン　右：f　をそれぞれ可視化したもの

   
   

反復法 - ヤコビ法 -

一般的なヤコビ法は、反復的に値を更新していくことで連立方程式を解く方法として知られています。

これから解説するヤコビ法は上記の一般的なヤコビ法とは異なり、時間方向と空間方向の双方に対して離散化した式を使って解いていきます。ですので、一般的なヤコビ法とは解いている式が異なりますので注意してください。

一般的なヤコビ法については話の本筋とそれてしまうためここでは解説しませんので、詳しく知りたい方は「ヤコビ法」で検索してみてください。

直接法でポアソン方程式を解くことができましたが、上記で分かるように直接法で解くためには巨大な行列を構築する必要があります。

2次元の 64×64グリッドの場合は (64*64)² = 16,777,216‬、
3次元の 64×64×64グリッドの場合は (64*64*64)² = 68,719,476,736

と途方もない行列の要素数になってしまいます。
そこで反復法というもので、省メモリでポアソン方程式を解いていきます。

ポアソン方程式を解く場合のヤコビ法では、以下の２次元ポアソン方程式の陽的差分の式（２次元ラプラス方程式の差分の式に f の項が追加されたもの）に対して、CFL条件である k ≦ h²/4 を課して時間方向に値を更新することで解いていきます。

\( \begin{align} \frac{\phi_{t+1}(x,y)-\phi_{t}(x,y)}{k} &= \frac{\phi_{t}(x+1,y)+\phi_{t}(x-1,y)+\phi_{t}(x,y+1)+\phi_{t}(x,y-1)-4\phi_{t}(x,y)}{h^2}-f(x,y) \\ \phi_{t+1}(x,y) &= \frac{\phi_{t}(x+1,y)+\phi_{t}(x-1,y)+\phi_{t}(x,y+1)+\phi_{t}(x,y-1)}{4}-\frac{h^2 f(x,y)}{4} \\   \end{align} \)

また、ポアソン方程式を反復法で解く場合の境界条件は、値が存在しない箇所（ボクセル）に任意の値を割り当てることで設定できます。

値 0 のデイリクレ条件（固定境界）の場合は、値が存在しない箇所（ボクセル）の値を 0 とすることで実装することができます。

Houdiniで実装

上記をHoudiniで実装していきます。

多くの並列での反復計算が必要なため OpenCL SOP でGPUを使い計算します。

1. OpenCL SOP で結果を書き込むための2Dボリュームを作成します。
   phi と __phi の2つのボリュームを作成するのは複数スレッドからの値の書き込みによる競合を回避するためです。
   

   

2. OpenCL SOP でヤコビ法を使い Aφ = f を解きます。
   メインカーネルで __phi に値を書き込み Write Back カーネルで __phi から phi に値をコピーすることで競合を回避します。反復回数は9000回です。
   

   

   

   ヤコビ法 ソースを展開 ソースを折りたたむ

   ---

   kernel void kernelName(
                    int phi_stride_x,
                    int phi_stride_y,
                    int phi_stride_z,
                    int phi_stride_offset,
                    int phi_res_x,
                    int phi_res_y,
                    int phi_res_z,
                    global float * phi ,
                    global float * __phi ,
                    global float * f
   )
   {
       int gidx = get_global_id(0);
       int gidy = get_global_id(1);
       int gidz = get_global_id(2);
       int idx = phi_stride_offset
               + phi_stride_x * gidx
               + phi_stride_y * gidy
               + phi_stride_z * gidz;
                
       float v1 = phi[phi_stride_offset
                    + phi_stride_x * max(gidx-1, 0)
                    + phi_stride_y * gidy
                    + phi_stride_z * gidz];
       v1 *= gidx != 0;
        
       float v2 = phi[phi_stride_offset
                    + phi_stride_x * min(gidx+1, phi_res_x-1)
                    + phi_stride_y * gidy
                    + phi_stride_z * gidz];
       v2 *= gidx != phi_res_x-1;
        
       float v3 = phi[phi_stride_offset
                    + phi_stride_x * gidx
                    + phi_stride_y * max(gidy-1, 0)
                    + phi_stride_z * gidz];
       v3 *= gidy != 0;
        
       float v4 = phi[phi_stride_offset
                    + phi_stride_x * gidx
                    + phi_stride_y * min(gidy+1, phi_res_y-1)
                    + phi_stride_z * gidz];
       v4 *= gidy != phi_res_y-1;
        
       float f_   = f[idx];
       float phi_ = (v1 + v2 + v3 + v4 - f_) / 4.0f;
       __phi[idx] = phi_;
   }
    
   kernel void writeBack(
                    int phi_stride_x,
                    int phi_stride_y,
                    int phi_stride_z,
                    int phi_stride_offset,
                    int phi_res_x,
                    int phi_res_y,
                    int phi_res_z,
                    global float * phi ,
                    global float * __phi ,
                    global float * f
   )
   {
       int gidx = get_global_id(0);
       int gidy = get_global_id(1);
       int gidz = get_global_id(2);
       int idx = phi_stride_offset
               + phi_stride_x * gidx
               + phi_stride_y * gidy
               + phi_stride_z * gidz;
                
       phi[idx] = __phi[idx];
   }
   

   ---

3. 結果を可視化すると問題なく解けていることが確認できます。
   

   

4. 直接法との誤差も ±1e-05 くらいの範囲に収まっていることが確認できます。
   
   

   左：直接法　右：ヤコビ法

反復法 - ガウス＝ザイデル法・SOR法 -

次に同じ反復法であるガウス＝ザイデル法とSOR法について実装していきます。
ここでも、一般的なガウス＝ザイデル法とSOR法とは解いている式が異なりますので注意してください。

ヤコビ法とガウス＝ザイデル法の違いは値の更新の仕方です。ヤコビ法では上下左右のセルの値は１つ前の反復の値を参照していましたが、ガウス＝ザイデル法では上下左右セルのうち２つのセルは現在の反復で更新した値を使用することで収束を早めます。

逐次計算の場合は全グリッドのループ中に φ の配列を上書きしていくだけで簡単に実装できますが、並列計算の場合は少し工夫が必要となります。並列計算ではチェス盤の赤マス・黒マスのように値を更新していくことでガウス＝ザイデル法を実現します。それにちなんでこの計算方法はRed-Black ガウス＝ザイデル法などと呼ばれます。

SOR法はガウス＝ザイデル法に加速パラメータ ω を追加することで、より収束を高速化させます。ω の値が1の場合はガウス＝ザイデル法と等価になります。

\( \phi_{i+1}(x,y) = \phi_{i}(x,y) + \omega(\tilde{\phi}_{i+1}(x,y) - \phi_{i}(x,y)) \\ \)

ということで、ガウス＝ザイデル法とSOR法をHoudiniで実装していきます。

Houdiniで実装

1. OpenCL SOPでガウス＝ザイデル法・SOR法で Aφ = f を解きます。
   メインカーネルで偶数セル、Write Back カーネルで奇数セルを更新していきます。
   

   

   ガウス＝ザイデル法・SOR法 ソースを展開 ソースを折りたたむ

   ---

   kernel void kernelName(
                    int phi_stride_x,
                    int phi_stride_y,
                    int phi_stride_z,
                    int phi_stride_offset,
                    int phi_res_x,
                    int phi_res_y,
                    int phi_res_z,
                    global float * phi ,
                    global float * __phi ,
                    global float * f ,
                    float  omega ,
                    float  invomega
   )
   {
       int gidx = get_global_id(0);
       int gidy = get_global_id(1);
       int gidz = get_global_id(2);
       int idx = phi_stride_offset
               + phi_stride_x * gidx
               + phi_stride_y * gidy
               + phi_stride_z * gidz;
                
       int solve = (gidx + gidy + gidz) % 2 == 0;
        
       float v1 = phi[phi_stride_offset
                    + phi_stride_x * max(gidx-1, 0)
                    + phi_stride_y * gidy
                    + phi_stride_z * gidz];
       v1 *= gidx != 0;
        
       float v2 = phi[phi_stride_offset
                    + phi_stride_x * min(gidx+1, phi_res_x-1)
                    + phi_stride_y * gidy
                    + phi_stride_z * gidz];
       v2 *= gidx != phi_res_x-1;
        
       float v3 = phi[phi_stride_offset
                    + phi_stride_x * gidx
                    + phi_stride_y * max(gidy-1, 0)
                    + phi_stride_z * gidz];
       v3 *= gidy != 0;
        
       float v4 = phi[phi_stride_offset
                    + phi_stride_x * gidx
                    + phi_stride_y * min(gidy+1, phi_res_y-1)
                    + phi_stride_z * gidz];
       v4 *= gidy != phi_res_y-1;
        
       float f_   = f[idx];
       float phi_ = phi[idx] * (!solve) + (invomega * phi[idx] + omega * (v1 + v2 + v3 + v4 - f_) / 4.0f) * solve;
       __phi[idx] = phi_;
   }
    
   kernel void writeBack(
                    int phi_stride_x,
                    int phi_stride_y,
                    int phi_stride_z,
                    int phi_stride_offset,
                    int phi_res_x,
                    int phi_res_y,
                    int phi_res_z,
                    global float * phi ,
                    global float * __phi ,
                    global float * f ,
                    float  omega ,
                    float  invomega
   )
   {
       int gidx = get_global_id(0);
       int gidy = get_global_id(1);
       int gidz = get_global_id(2);
       int idx = phi_stride_offset
               + phi_stride_x * gidx
               + phi_stride_y * gidy
               + phi_stride_z * gidz;
                
       int solve = (gidx + gidy + gidz) % 2 != 0;
        
       float v1 = __phi[phi_stride_offset
                      + phi_stride_x * max(gidx-1, 0)
                      + phi_stride_y * gidy
                      + phi_stride_z * gidz];
       v1 *= gidx != 0;
        
       float v2 = __phi[phi_stride_offset
                      + phi_stride_x * min(gidx+1, phi_res_x-1)
                      + phi_stride_y * gidy
                      + phi_stride_z * gidz];
       v2 *= gidx != phi_res_x-1;
        
       float v3 = __phi[phi_stride_offset
                      + phi_stride_x * gidx
                      + phi_stride_y * max(gidy-1, 0)
                      + phi_stride_z * gidz];
       v3 *= gidy != 0;
        
       float v4 = __phi[phi_stride_offset
                      + phi_stride_x * gidx
                      + phi_stride_y * min(gidy+1, phi_res_y-1)
                      + phi_stride_z * gidz];
       v4 *= gidy != phi_res_y-1;
        
       float f_   = f[idx];
       float phi_ = __phi[idx] * (!solve) + (invomega * __phi[idx] + omega * (v1 + v2 + v3 + v4 - f_) / 4.0f) * solve;
       phi[idx] = phi_;
   }
   

   ---

   
   
2. 結果を可視化すると問題なく解けていることが確認できます。
   

   

3. SOR法による収束は高速で、ヤコビ法では反復回数 9000回かかっていた計算が、ω の値が 1.9 の場合 反復回数 200回で直接法との誤差が ±1e-05 くらいの範囲に収まります。
   

   
   

グリーン関数

ポアソン方程式を解くのにグリーン関数（高速フーリエ変換）を使う方法もあります。今回は、重力ポテンシャルを求めてみます。

以下の手順で重力ポテンシャルを求められます。

3次元重力場のポアソン方程式は以下です。（G は重力定数、ρ は密度分布）

\( \Delta\Phi = -4{\pi}G{\rho} \\ \)

3次元ポアソン方程式のグリーン関数の定義は以下です。

\( G(r,r^{\prime}) = \frac{1}{4{\pi}|r-r^{\prime}|} \\ \)

重力ポテンシャル Φ は、密度分布 ρ とグリーン関数 G(r) を畳み込み積分することで求められます。

\( \begin{align} \Phi(r) &= -4{\pi}G\int\rho(r^{\prime})\frac{1}{4{\pi}|r-r^{\prime}|}dr^{\prime} \\ &= -4{\pi}G(\rho*G)(r) \\ \end{align} \)

畳み込み積分はフーリエ変換の結果を乗算することで求められます。F はフーリエ変換です。

\( \mathscr{F}(\rho*G) = \mathscr{F}(\rho)\mathscr{F}(G) \\ \)

よって、重力ポテンシャルは 密度分布 ρ をフーリエ変換した値 ρ~ と、グリーン関数 G をフーリエ変換した値 G~ と、-4πG を掛け合わせた値（Φ~）を逆フーリエ変換することで求められます。

\( \tilde{\Phi}(k) = -4{\pi}G\tilde{\rho}(k)\tilde{G}(k) \\ \)

また、上記の方法で何もせずにポアソン方程式を解いた場合の境界条件は周期境界条件となります。

今回は、孤立境界条件で解くためにゼロパディングを使用します。

下の図のように3Dボリュームを8分割したうちの1つの象限に密度分布を配置し、他の象限の密度分布を0とするとこで孤立境界条件で解くことができます。

Houdiniで実装

上記をHoudiniで実装していきます。

1. 3Dボリュームを作成します。
   高速フーリエ変換を行うため、実部と虚部の値を格納するボリュームが2つ必要となります。
   

   

2. グリーン関数を作成します。
   

   

3. 密度分布を作成します。
   

   

4. グリーン関数と密度分布を順方向3次元高速フーリエ変換します。
   

   

5. Φ~（Φチルダ）を計算します。
   グリーン関数は遇関数のため、グリーン関数のフーリエ変換の結果の虚部の値（奇関数の成分）は誤差として無視してしまします。
   

   

6. Φ~（Φチルダ）を逆方向3次元高速フーリエ変換を行い Φ を求めます。
   

   

7. 結果を可視化します。
   

   

8. 重力ポテンシャルを総和で求めた結果と比較します。
   計算対象の象限に関して、問題なく解けていることが確認できます。
   
   

   左：総和　右：グリーン関数

まとめ

今回はHoudiniでポアソン方程式を色々な方法で解いてみました。

Houdiniでは並列計算が容易に実装できることもあり、並列化できる方法であれば、簡単に高速で動作する処理を作成することができます。

ノード単位で処理を分割することで、処理の順番を簡単に入れ替えたり、処理を行った結果を簡単に視覚化することもできます。

ですので、著者はHoudiniを数式や論文を実装するにあたってのプロトタイピングツールとしても使ったりしています。

また、今回ポアソン方程式を解くいろいろな方法を紹介しましたが、方法による優劣はなく目的によって最適な方法を使い分けるのが重要だと著者は考えています。

例えば、リアルタイムでのシミュレーションが目的で解の精度を必要としないなら、ヤコビ法などで数回の反復で計算を打ち切ってしまうのも手です。

正確な解が欲しくてメモリに余裕があるようでしたら直接法で解くのもまた手です。

マシンのメモリや、計算速度、計算精度などと相談して、ケースバイケースで最適な方法を探ってみてください。

この記事がポアソン方程式を解く機会があったときの手助けになれば幸いです。