NF前処理編をお送りします。極く一般的なILU前処理に比べ格段に優れた発想とつくづく思う。俺的には最も気に入っている前処理法。考案者に乾杯です。This is the practical lecture course about the NF preconditioning method. A comparison with a very common ILU preconditioning would lead me to reassure the NF preconditioning extremely superior method. This is my favorite preconditioning. I wish to toast the the inventor for creating this superb factorization algorithm.

Overall NF Algorithm

下図４.０に示す一般的な三次元問題\((x,y,z)=(4,3,2)\)に於ける係数行列は\(A\)は次式で示される。

$$A = D + L_{1}+ U_{1} + L_{2} + U_{2} + L_{3} + U_{3}$$

ここで、\(D\)は対角成分、\(U\)は上三角対角行列、\(L\)は下三角対角行列である。又、添字は対応する流れの次元を示す。

今､係数行列\(A\)をブロック三重対角構造(Block Tridiagonal Structure)に分割してみよう。係数行列\(A\)において､最も外側にある帯行列\(L_{3}\)と\(U_{3}\)は面(Plane)間の繋がり(５･６点有限差分展開)を示しており、面毎に面ブロック三重対角構造に分割できる｡更に、各々の面は線(Line)ブロック三重対角構造に再分割できる｡その再分割されたブロック三重対角構造は線間の繋がり(３･４点有限差分展開)を示す帯行列\(l_{2}\)と\(U_{2}\)から構成される｡又､線ブロックは各々\(l_{1}\)と\(U_{1}\)からなる小単位の対角構造からなる。このようにブロック三重対角構造は下成分に上行列から帯行列を一つ導入し､上成分には下行列から帯行列を一つ導入する事によって形成されうる。

依って､係数行列\(A\)は次のような三次元ブロック三重対角行列(Block Tridiagonal Matrix)\(M\)に近似できる。二次元面(Plane)結合で構成されるブロック三重対角行列\(P\)、一次元ブロック三重対角行列\(T\)、如いては対角行列\(G\)にて、再帰的に分解近似できる。

$$\begin{align}
M &= (P+L_{3})(I+P^{-1}U_{3}) \\
P &= (T+L_{2})(I+T^{-1}U_{2}) \tag{4.0} \\
T &= (G+L_{1})(I+G^{-1}U_{1}) \\
G &= D – L_{1}S^{-1}U_{1} – colsum(L_{2}T^{-1}U_{2}) – colsum(L_{3}P^{-1}U_{3})
\end{align}$$

ここで、\(I\)は単位行列(Identity Matrix)である。

対角行列\(G\)は\(colsum(L_{2}T^{-1}U_{2}+L_{3}P^{-1}U_{3})\)を差引く事によって、誤差行列(Error Matrix)の列(column)の総和は零となる。ここで、\(colsum\)とはcolumn毎に行列成分の総和をとったものである。詰まり、元の係数行列\(A\)に対する誤差行列\((M-A)\)は、\(L_{2}T^{-1}U_{2} – colsum(L_{2}T^{-1}U_{2}) + L_{3}P^{-1}U_{3} – colsum(L_{3}P^{-1}U_{3})\)で与えられ、Columnの総和\((colsum)\)を零にする設定であることが分かる。

以上の近似式を使って､連立1次方程式\(Mx=r\)を次のように再帰的に分解していく。

具体的には、三重対角行列\(M\)は次の２組の方程式を解く事によって行う。

$$\tag{4.1}(P+L_{3})(I+P^{-1}U_{3}) x=r$$

$$\left\{
\begin{align}
(P+L_{3})w &=r \\
(I+P^{-1}U_{3})x &= w
\end{align}
\right.$$

本例では、\((P+L_{3})\)行列と\((I+P^{-1}U_{3})\)行列は下図４.１に示す２x２行列である。

最初の二次元面構造から最後の二次元面構造に亘って、第１式によって\(w=P^{-1}(r-L_{3}w)\)から前進消去にて\(w\)を求める。次に、得られた\(w\)をもとに、最後から最初の二次元面構造まで遡り、第２式によって\(x=w-P^{-1}U_{3}x\)から後進代入法にて\(x\)を求める。

さて､上記過程において\(s=P^{-1}q\)を算定する必要がある。\(P\)は次式で与えられる。

$$\tag{4.2}(T+L_{2})(I+T^{-1}U_{2})s =q$$

$$\left\{
\begin{align}
(T+L_{2})w &=q \\
(I+T^{-1}U_{2})s &= w
\end{align}
\right.$$

ここで、 \((T+L_{2})\)行列と\((I+T^{-1}U_{2})\)行列は下図４.２で示す３x３行列である。

先と同様に、面内の一次元線結合構造から最後の一次元線結合構造に亘って、第１式\(w=T^{-1}(q-L_{2}w)\)から前進消去にて(w)を求める。次に、得られた\(w\)をもとに、最後から最初の一次元線結合構造まで遡って、第２式によって\(s=w-T^{-1}U_{2}s\)から後進代入法にて\(s\)を求める。本例に於ける\(T_{a},T_{b},T_{c}\)は最小の一次元ブロック単位(４x４行列)である。

ここで、上記過程において\(v=T^{-1}z\)を算定しなければならないが、\(T\)は次式で与えられる。

$$\tag{4.3}(G+L_{1})(I+G^{-1}U_{1})v =z$$

$$\left\{
\begin{align}
(G+L_{1})w &=z \\
(I+G^{-1}U_{1})v&= w
\end{align}
\right.$$

先と同様の過程によって、最終的に一次元\(T^{-1}z\)を第１式を前進消去法、第２式を後進代入法によって其々解くのだが、これは一般的な対角行列を解く事に他ならない。

解の初期値\(X^{0}\)は次式で与える。

$$X^{0} = M^{-1}b \tag{4.4}$$

以上､NF前処理はブロック三重対角構造の係数行列に於いて外側から内側の帯行列へ向かって順次再帰分解していく手法である。ILU前処理に比べて、NF前処理はStraight Forwardである。各\(Colsum\)の算定、各\(M,P,T\)三重対角行列の解法を逐次実行することで求解できる。

Solution of Diagonal Matrix \(G\)

対角行列\(G\)の算定方法について説明する｡

先のプロセスとは全く逆のプロセスを辿って解く事になる｡
詰まり、格子単位から線､面構造へと順次外側の三重対角行列構造へ逆算していく。

手順は以下。

１） 格子毎に\(Ll_{1}G^{-1}U_{1}\)を順次求め､線ブロック構造内の\(L_{1}G^{-1}U_{1}\)を求める。

２）\(colsum(L_{2}T^{-1}U_{2})\)の算定

　　(1)\(L_{2}T^{-1}\)は次式で与えられる。
$$x=L_{2}T^{-1} \rightarrow Tx=L_{2}$$

　　(2)\(T\)は次の三重対角行列を解くことで容易に求まる。
$$(G+L_{1})(I+G^{-1}U_{1})x=L_{2}$$

　　(3)\(x\)が求まれば、\(L_{2}T^{-1}U_{2}=U_{2}x\)より、\(L_{2}T^{-1}U_{2}\)が求まる。

　　(4)\(L_{2}T^{-1}U_{2}\)の転置に単位行列を掛合わせ、\(colsum(L_{2}T^{-1}U_{2})\)を求める。

３） \(colsum(L_{3}P^{-1}U_{3})\)の算定

　　上記\(T\)と同様の手法で求める。

以上のプロセスを全ての面構造に対して繰り返し、対角行列\(G\)を求める。

Closing

NF前処理の収束性は軸方向の順序付けによって大きく左右される事が知られている。
特に、Transmissibility(易動度)の影響が重大で、易動度が大きい軸ほど内側の帯行列になるように軸方向を定めるべきとの研究結果がある。従って、順序付けに際しては、易動度、面構造及び線構造体数を総合的に留意しなければならない｡

実際、NF前処理の計算過程に於けるMaterial Balance ErrorはILU処理に比べ優れている。しかしながら、”ECLIPSE”に於けるNewton-Raphson法の収束判定基準は”IMEX”に比べて厳しい基準が課されていて、”IMEX”基準に緩和させることで一層の収束性の向上を図ることが可能である。

次回のテーマは”SMARTSOL”の特記事項の順応陰解法(AIM)及び強結合坑井モデル(Strongly Coupled Well Model)に関する座学主体の講義を予定している。

以上
つづく